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Chapter  1 
Introduction 


1.1  Purpose  of  this  report 

This  report  is  focused  on  the  numerical  simulation  of  the  micro-mechanical  response  of  concrete  under 
dynamic  tensile  loading.  The  microstructure  plays  an  important  role  in  the  behavior  of  heteogeneous 
materials  under  impact  loading.  For  instance,  concrete  has  a  higher  resistance  under  impact  loading  than 
under  static  loading.  This  is  partly  due  to  the  microstructure.  The  micro-mechanical  background  of  this 
observation  is  the  main  topic  of  this  report. 


1.2  Relevant  concepts 

1.2.1  Representative  volume  element 

In  this  investigation,  a  microstructure  or  Representative  Volume  Element  (RVE)  is  used.  An  RVE  for  a 
material  point  of  a  continuum  is  a  material  volume  which  is  sufficiently  representative  of  the  infinitesimal 
material  neighbourhood  of  that  material  point  (cf.  [13]).  It  is  regarded  as  a  heterogeneous  material  with 
spatially  varying,  but  known,  constitutive  properties.  It  is  of  importance  that  statistical  properties  of 
the  state  variables  are  considered  invariant  of  the  position  in  the  material.  This  is  called  statistical 
homogeneity. 

Since  concrete  is  a  highly  heterogeneous  material,  the  size  and  shape  of  the  RVE  should  be  chosen 
properly  with  respect  to  the  dimensions  of  concrete  structures.  For  example,  a  minimum  aspect  ratio, 
which  is  the  ratio  between  the  microstructure  dimension  and  the  maximum  grain  size,  should  be  taken 
of  3  to  5  (cf.  [14]).  Another  essential  step  in  numerical  homogenization  is  the  appropriate  choice  of  the 
microscopic  boundary  condition.  Since  the  homogenization  procedure  is  not  followed,  this  issue  is  not 
addressed  in  this  investigation. 


1.2.2  Softening  and  localization 

The  main  cause  of  the  nonlinearity  of  concrete  is  cracking,  which  is  primarily  due  to  the  limited  capacity 
of  concrete  to  sustain  tensile  stresses  (  or  perhaps  better,  tensile  strains).  Three  main  approaches  can 
be  distinguished  for  the  modelling  of  cracking  in  concrete,  namely  discrete  crack  models,  for  example 
the  cohesive  surface  methodology  adopted  by  Tijssens  (cf.  [12])  for  simulation  of  the  fracture  process  of 
concrete  at  macro-  and  meso-level,  the  smeared  crack  models  and  lattice  models.  Within  the  smeared 
crack  concept,  a  cracked  solid  is  imagined  to  be  a  continuum  where  the  notions  of  stress  and  strain  remain 
valid.  Hence,  continuum  constitutive  models  can  be  used  (cf.  [3]). 

Concrete  can  be  classified  as  a  softening  material  which  shows  a  reduction  of  the  load  carrying  capacity 
accompanied  by  localized  deformations  after  reaching  the  limit  load.  Experimentally  obtained  load- 
displacement  characteristic  exhibits  a  descending  branch  which  is  commonly  known  as  structural  softening 
behavior.  Structural  softening  as  global  behavior  is  due  to  micro-cracking  in  concrete.  At  a  certain  stage 
of  growth  and  coalescence  of  microcracks,  deformations  grow  in  a  small  portions  of  the  material  and  the 
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rest  of  the  structure  tends  to  unload.  This  phenomenon  is  known  as  strain  localization  and  eventually 
leads  to  the  formation  of  macroscopic  cracks. 

For  a  continuum  description  of  structural  softening  behavior  or  strain  localization  phenomenon,  a 
material  law  with  a  negative  slope  after  the  limit  load  is  used.  This  is  called  a  strain  softening  material 
law. 

1.2.3  Regulairization  methods 

Unfortunately  the  field  equations  that  describe  the  motion  of  the  body  lose  hyperbolicity  and  become 
elliptic  as  soon  as  strain  softening  occurs,  when  a  standard  rate  independent  continuum  model  is  combined 
with  a  strain  softening  law.  The  domain  is  split  up  into  an  elliptic  part,  in  which  waves  have  imaginary 
wave  speeds  and  are  not  able  to  propagate  {standing  waves),  and  into  a  hyperbolic  part  with  propagating 
waves.  The  initial  value  problem  becomes  ill-posed  and  can  no  longer  be  a  proper  description  of  the 
underlying  physical  proplem.  Localization  zones  stay  confined  to  a  line  with  zero  thickness  (or  a  discrete 
plane  in  a  three  dimensional  continuum).  The  energy  consumed  in  the  failure  zones  is  zero.  These  two 
results  are  in  contradiction  with  experimental  results.  In  numerical  simulation,  the  finite  element  tries 
to  capture  the  localization  zone  with  zero  thickness  which  results  in  a  mesh  sensitivity.  The  reason 
for  this  fact  is  that  the  adopted  softening  laws  are  deduced  phenomenologically  and  contain  no  actual 
information  on  the  microscopic  behavior  of  the  material,  which  plays  an  important  role  during  failure 
process  (cf.  [11]). 

In  order  to  properly  model  the  softening  solid  within  continuum  framework,  a  few  approaches  have  been 
proposed  to  regularize  the  governing  equations,  such  as  fracture  energy  based  models,  rate  dependent 
models,  nonlocal  models,  gradient  models,  micro-polar  models  (e.g.  Cosserat  model),  and  embedded 
discontinuity  models.  The  common  feasure  of  these  models  is  the  inclusion,  implicitly  or  explicitly,  of 
a  length  scale  into  the  field  equations  or  the  description  of  the  material.  With  these  approaches,  well- 
posedness  of  the  initial  value  problem  can  be  preserved,  the  governing  equations  do  not  necessarily  lose 
hyperbolicity,  and  a  solution  with  real  wave  speeds  can  be  obtained  at  the  onset  of  strain  softening.  The 
numerical  results  converge  to  a  finite  size  of  the  localization  zone  and  nonzero  energy  consumption  in  the 
failure  process  zones  upon  mesh  refinement.  The  mesh  dependent  solution  can  be  avoided. 

In  this  report,  the  implicit  gradient  dependent  damage  model  and  two  rate  dependent  viscoplastic 
models  will  be  adopted  as  the  constitutive  descripion  of  the  RVE.  A  length  scale  is  incorporated  in  these 
models  which  preserve  mesh  independence. 


1.3  Arrangement  of  the  contents 

Chapter  2  starts  with  basic  equilibrium,  constitutive  and  kinematic  equations  for  the  motion  of  an  inelastic 
body.  The  finite  element  discretization  of  the  virtual  work  equation  and  implicit  time  integration  scheme 
will  be  derived. 

Chapter  3  introduces  two  kinds  of  continuum  models  used  in  this  report,  the  continuum  plasticity 
model  and  the  continuum  isotropic  damage  model.  The  standard  elastoplastic  model,  Perzyna  viscoplastic 
model,  Duvaut-Lions  viscoplastic  model,  the  implicit  constant  gradient  enhanced  damage  model  and  the 
strain-based  transient-gradient  damage  model  are  described  including  the  algorithmic  treatment. 

In  Chapter  4  the  numerical  results  of  a  one-dimensional  simulation  with  the  models  introduced  in 
Chapter  3  are  given  and  the  regularization  effects  are  addressed. 

In  Chapter  5  the  results  of  the  numerical  simulation  of  the  dynamic  response  of  concrete  at  meso-level 
are  presented.  The  characteristics  of  each  model  are  illustrated  and  compared. 

The  report  is  finished  with  conclusions  and  a  discussion  in  Chapter  6. 


Chapter  2 

Formulation  of  the  initial  value  problem 


The  basic  structure  of  the  initial  value  problem  (IVP)  for  nonlinear  solids  is  outlined  with  a  restriction 
to  small  deformation  gradients.  We  use  a  finite  element  representation  of  the  virtual  work  equation  to 
derive  the  semi-discrete  initial  value  problem.  A  full  discretization  is  obtained  if  the  time  integration 
has  been  carried  out.  An  incemental-iterative  Newton-Raphson  procedure  is  used  to  solve  the  resulted 
nonlinear  equations.  Throughout  this  report  matrix-vector  notation  is  used. 


2.1  Preliminary  equations 


We  consider  a  body  B  with  volume  and  surface  dU  —  dClt  where  dQt  and  dQu  are  traction  and 

kinematic  boundaries,  respectively,  with  dClt  ^  =  0.  The  equations  of  motion  without  damping,  the 

kinematic  equations  and  the  corresponding  initial  and  boundary  conditions  read 


L^cr  -1-  p  =  Rii 

in 

n 

(2.1) 

e  =  Lu 

in 

n 

(2.2) 

u(x,  0)  =  uo 

v(x,  0)  =  vo 

in 

n 

(2.3) 

u(x,  t)  =  u 

on 

dClu 

(2.4) 

t  =  N^a 

on 

dQt 

(2.5) 

where  the  differential  operator  L  is  defined  as 


( 

0 

^  0  0 

0  1  0 

0  lx  1 

1  0  ly 

dz  ' 
0 

1  j 

(2.6) 

Stress  and  strain  vectors  cr  and  e  are 

represented  as 

II 

7xx  ^yy 

^zz  CTxy 

<7yz 

^zx\ 

(2.7) 

6^  =  [e 

XX  ^yy 

€.zz 

2eyx 

2^2  x] 

(2.8) 

u,  V,  p,  t  represent  the  displacement, 
the  outward  normal  n  defined  as 

velocity,  body  force  and  surface  traction  with  matrix  N  related  to 

!  Tlx  0 

0  riy 

\  0  0 

0  %  0 

0  Tlx  Tlx 

Tlx  0  Tly 

Tlx  \ 

Tlx  J 

(2.9) 

3 


4 


The  density  matrix  R  is  equal  to  diag[/9  p  p]  with  p  the  material  density. 

The  constitutive  euqations  can  be  given  in  the  general  rate  format 

&  =  (2.10) 

with  matrix  0“=  containing  the  tangent  stiffness  moduli.  Two  types  of  material  constitutive  descriptions 
are  used  in  this  report,  the  continuum  damage  model  and  the  continuum  plastic  model.  They  will  be 
detailed  in  Chapter  3. 


2.2  Weak  form  of  the  initial  value  problem 


While  Eqs.  2.1  to  2.5  describe  the  motion  of  a  body  in  a  strong  sense,  a  weak  form  of  these  equations  is 
obtained  by  setting 


f  (5u^[L^<7  +  p  —  RiiJdQ  =  0 

Jn 


(2.11) 


in  which  <5  denotes  the  variation  of  a  quantity.  With  the  aid  of  Green’s  theorem,  this  equation  can  be 
transformed  into  (cf.  [11,  15,  3,  1]) 


f  5u^[Ru]dr2+  f  f  5u'^pdCl+  f  (Ju^tdF 

Jn  Ja  Ja  Jerit 


2.3  Finite  element  discretization  of  the  I  VP 


(2.12) 


For  the  use  of  implicit  time  integration  scheme,  Eq.  2.12  is  considered  to  be  valid  at  time  t+At.  The 
updated  stress  at  global  (I+l)-th  iteration  reads 


'^t+At  —  <^4  +  ~  ^t+At  + 

in  which  the  correction  of  the  (I)-th  incremental  stress  dcr^^^  can  be  written  as 


(/) 


(2.13) 


(2.14) 


t+A< 


Then  the  (I+l)-th  iteration  form  of  Eq.  2.12  can  be  expressed  as 
Jn  Jq  Jn 


I 

JdQ.t 


-  f 

Jq 


n 


(2.15) 


Making  use  of  the  standard  finite  element  formulation  the  continuous  displacement  field  u  and  the 
continuous  acceleration  field  u  are  discretized  as 


u(x,  t)  —  H(x)a(t) 
u(x,  t)  =  H(x)a(t) 
e(x,t)  =  B(x)a(t)  =  LHa(t) 

where  H(x),  B(x)  the  interpolation  matrix  and  the  strain-nodal  displacement  matrix. 
Substitution  of  Eqs.  2.16  to  2.18  and  2.14  into  Eq.  2.15  leads  to 

(5a^  f  B.'^RKa[^_^^^dQ  +  5a^  [  =  da^  [ 

Jq  Jq  Jq 

7 

JdQt 


(2.16) 

(2.17) 

(2.18) 


4-(5a^ 


H^tli+Vdr  -  5a'^  f 
Jn 


't+At^ 


(2.19) 
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Because  this  identity  must  hold  for  any  admissible  da  we  obtain 

(2.20) 

where  the  subscript  t+At  is  dropped  for  simplicity.  The  mass  matrix  M,  the  stiffness  matrix  the 
internal  force  vector  and  the  external  load  vector  given  by 

M  =  /  H^RHdn 

Jn  Jn 

c=  [  4^+^)=  f  f  (2.21) 

Jn  Jn  JsQt 

Eq.  2.20  represents  the  semi-discrete  nonlinear  equation  of  motion  governing  the  response  of  a  system  of 
finite  elements. 


2.4  Time  integration  scheme 


After  spatial  integration  of  the  virtual  work  equation  we  apply  the  direct  integration  method  to  obtain  a 
full  discrete  equation  of  motion.  A  one  step  time  integration  scheme  of  Newmark  family  is  used  in  this 
report,  in  which  quadratic  expansion  is  used 


St+At  =  at  -f  at  At  -I-  [(—  —  /?)at  -f-  /3at+At]At^ 

at+A<  =  at  -I-  [(1  —  7)at  4-  7at+At]At 

'1  1' 


(2.22) 


(2.23) 


We  use  (/3,7)  =  ( j)  for  all  calculations  performed  in  this  report,  which  is  called  the  average  acceleration 


method  or  trapezoidal  rule.  We  can  rewrite  Eq.  2.22  as 

^+aV  =  coAa^^+^^  -  ciat  -  C2at 
where  the  coefficients  are  given  by 


Co  = 


pAt^ 


Cl 


1 

jAt 


2p  ^ 


The  (I-t-l)-th  incremental  displacement  can  generally  be  expressed  as 
Aa(^+i)  =  Aa(')+da(') 

Then  we  obtain 

=  coda^-'^^  -f  coAa^^^  -  Cjat  -  cadt  =  coda^^^  4-  ^2At 
Substitution  of  Eq.  2.27  into  Eq.  2.20  leads  to  the  following  system  of  equations 
K(/)da(^)  = 

where  the  dynamic  tangent  matrix  and  the  effective  load  vector  are  defined  as 

K(^)  =  K<^)  4-  coM 

f(/+i)  f(/+i)  _ 


(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

(2.30) 


respectively. 

The  algorithm  for  implicit  integration  of  the  semi-discrete  equation  of  motion  is  outlined  in  Table 
2.1.  This  scheme  corresponds  to  a  full  Newton- Raphson  scheme  where  the  tangential  stress-strain  matrix 
is  updated  at  each  iteration. 


1.  calculate  constants 


^0  — 

*^1  ■“  PAt 

02—2/3  ^ 

03  =  (1  -  7)At 

<1 

f- 

II 

2.  initialize  a(0), 

a(0). 

a(0). 

t  =  0  s 

3.  initialize  Aa(°\ 

kw, 

K(°), 

f(i) 

ext  1 

f(0)  f(l) 

4.  solve  for  the  delta-incremental  displacements 


5.  update  the  iterative  incremental  displacements,  accelerations,  velocities  and  displacements 


^+At  > 


'‘t+At  ’ 


^t+At 


6.  check  convergence  criterion,  if  not  converged:  I  < —  I-l-l 

update  K(^),  KW,  Go  to  4. 


7.  set  new  accelerations,  velocities  and  displacements. 

r  .  ..  1  r  (7+1)  -  (7+1)  -(7  +  l)l 

[at+A«>  St+Af,  a^+AtJ  —  l^+At  ’  ^t+At  >  ^+At  J 


Table  2.1:  Global  implicit  integration  procedure 


Chapter  3 

Constitutive  models 


In  this  Chapter,  the  material  models  used  to  describe  the  mechanical  response  of  concrete  are  discussed. 
The  corresponding  algorithmic  treatments  are  also  given  for  implementation  in  a  finite  element  pro¬ 
gramme. 


3.1  Plasticity  model 

Plasticity  models  with  an  appropriate  yield  function  and  flow  rule  can  be  used  to  describe  mechanical 
behavior  of  concrete  under  compressive  loading  and  tensile  loading.  For  example,  the  yield  functions  of 
Von  Mises  and  Drucker-Prager  can  be  used  for  model-II  dominated  failure  problems.  On  the  other  hand 
the  Rankine  criterion  is  suited  to  predict  a  mode-I  failure  pattern. 

We  will  use  three  plasticity  models  in  the  simulations,  namely  the  standard  Von  Mises  plastic  model, 
Perzyna  viscoplastic  model  and  Duvaut-Lions  viscoplastic  model.  As  illustrated  in  Chapter  1,  the  rate 
dependence  of  plastic  flow  characterizing  the  latter  two  models  implicitly  introduces  an  internal  length 
scale  in  the  constitutive  descriptions  and  preserves  the  well-posedness  of  the  IVP  upon  the  onset  of  strain 
softening. 


3.1.1  Standard  plasticity  model 

In  the  flow  theory  of  plasticity,  the  material  behavior  can  be  written  in  a  general  rate  format 

k  =  i-e  +  €p  (^-l) 

d-  =  D"(e-€p)  =  D=ee  (3-2) 

kp  =  Am  (3-3) 

where  D®  is  the  elastic  stiffness  tensor,  A  is  the  plastic  multiplier  and  m  is  the  direction  of  plastic  flow. 
Loading  and  unloading  conditions  can  be  conveniently  expressed  by  Kuhn-Tucker  conditions 

/<0,  A>0,  /A  =  0  (3.4) 


where  f(<T,K)  is  the  yield  function  with  a  number  of  internal  variables  collected  in  a  vector  k  which 
describe  the  plastic  strain  history. 

For  isotropic  hardening/softening  this  vector  reduces  to  a  scalar-valued  hardening/softening  parameter 
K  depending  on  strain  history  through  invariant  measures  of  the  plastic  strain  €p.  According  to  the  strain 
hardening  hypothesis 

K=^2/3(ep)^Q€p  (3.5) 

where  Q  =  diag[l,l,l,l/2,l/2,l/2].  This  definition  of  k  is  used  in  this  report. 
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We  make  use  of  a  Von  Mises  yield  function  which  is  described  as 

/(cr,  k)  =  -  a-(/c)  (3.6) 

where  J2  is  the  second  invariant  of  the  deviatoric  stress.  The  combination  of  flow  theory  of  plasticity 
with  Von  Mises  yield  contour  is  also  named  the  J2-flow  theory. 

The  isotropic  softening  laws  have  been  utilized, 


-  the  linear  softening  law 


aiK)  = 


EnK  +  <70  0  <  K  <  Ku 

0  K>  K-a 


(3.7) 


where  jB„  is  the  softening  modulus,  do  the  initial  yield  stress  and  the  ultimate  value  of  n  at  which  d 
has  reduced  to  zero; 


-  the  exponential  softening  law 

a{K)  =  doo  +  (d-Q  -  doo)e^“hK)  (3  gj 

where  doo  is  the  yield  stress  at  infinite  value  of 
Another  condition  for  plastic  straining  is 


/  =  0 


(3.9) 


which  is  usually  called  Prager’s  consistency  condition. 

On  the  material  (local)  level,  the  constitutive  rate  equation  should  be  integrated  to  calculate  the  stress 
o’t+At  at  time  t+At  and  to  provide  the  tangential  material  stiffness  matrix  (cf.  [3]).  This  is  achieved  by 
utilizing  the  fully  implicit  Euler  Backward  integration  scheme,  i.e.  the  closest  point  projection  algorithm. 
The  incremental  form  of  constitutive  equation  Eq.  3.1  to  Eq.  3.3  is 


t+At  =  CTt  +  Act 

(3.10) 

Act  =  D®(Ae  —  Acp) 

(3.11) 

Aep  =  AAmf+A« 

(3.12) 

We  cast  Eq.  3.10  to  Eq.  3.12  in  the  form 

g(o't+At>  At+At)  —  0 

(3.13) 

where 

g(o’t+At)  At+A«)  —  ^t+At  —  'D^f-t+At  +  AAD'^nit+At 

(3.14) 

To  determine  (Xt+At  and  AA  we  must  augment  Eqs.  3.13  by  the  requirement  that  the  yield  condition  is 
complied  with  at  the  end  of  a  loading  step 

fi<Xt+At,  X't+At)  =  /(o’ t+At,  At+At)  =  0  (3.15) 

A  Newton-Raphson  procedure  is  set  up  at  integration  point  level  to  solve  this  nonlinear  system  of 
equations 


- 1 

+  + 

_ 1 

> 

_ 1 

■  is. 

da- 

1 

a\ 

-1  p 

.  ^ttL  . 

< 

at 

dX  J 

S{cr'^+Av>^t+At)  ' 
fi^t+At^^t+At)  . 


(3.16) 
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where  k  is  the  iteration  counter  of  the  local  Newton- Raphson  iteration.  In  view  of  Eq.  3.14  the  differentials 
in  Eq.  3.16  can  be  elaborated  as 


|S  =  a  =  i+aai>-^ 

da  da 

%  T-xP  -  -rxp/  ..3k  5m, 

^  =  D^m  =:  D^(m  +  AA--)  =  D^(m  +  AA-— ) 


(3.17) 

(3.18) 


A  consistent  linearization  of  Eq.  3.10  to  Eq.  3.12  together  with  consistency  condition  Eq.  3.9  leads  to 
the  algorithmic  tangential  stiffness  relation 


D^Pe 


U- 


Uffm^U 
h  +  n^Uin 


(3.19) 


where  U  =  A  h  =  =  §§■ 


3.1.2  Perzyna  viscoplastic  model 

In  the  viscoplastic  model  according  to  Perzyna  ([10])  the  viscoplastic  strain  rate  is  defined  similar  to  the 
plastic  strain  rate  in  standard  plasticity  (i.e.  Eq.  3.1  and  Eq.  3.2)  with  a  modified  definition  of  the  plastic 
multiplier. 


e  =  Am 
A  =  77<A(/) 


(3.20) 

(3.21) 


in  which  77  is  the  viscosity  parameter,  (j)  is  an  arbitrary  function  of  the  yield  function  for  which  we  take 
a  power  law  expression  according  to  (cf.  [11,  15]) 


(3.22) 


l\^o/ s 

with  N  a  constant,  do  the  initial  yield  stress  and  = 

A  full  implicit  integration  scheme  is  used  to  obtain  the  stress  Cj+At  at  time  t-fAt  and  the  consistent 
tangential  stiffness  matrix  The  same  incremental  form  of  constitutive  equations  as  Eq.  3.10  to 

Eq.  3.12  can  be  obtained,  then  we  obtain  the  following  equations 


2 


AA 

r  =  4>{(Jt+At,  At+At)  —  =  0 

Aa  =  D®[Ae  -  AAmt+At] 


(3.23) 

(3.24) 


A  local  Newton-Raphson  procedure  is  applied  to  compute  AA  and  Aa.  During  the  linearization  a 
consistent  tangent  stiffness  matrix  is  obtained 


=  P  -  -P 
a 


m  -t-  AA 


3m 


d(f) 


da 


in  which 


(D®)-^-h  AA 


3m 

da 


1-1 


a  = 


'  d4>' 

T 

P 

A  X 

m-l- AA-— 
3A 

-1- 


1  _  ^ 
■qAt  3A 


(3.25) 

(3.26) 

(3.27) 


Algorithmic  treatment  for  Perzyna  viscoplastic  model  is  outlined  in  Table  3.1. 
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1.  update  strain  increment  and  total  strain  Ac,  e 

2.  compute  trial  stress  cr*’’  =  D®(e  -  c^p) 

3.  if  >  0  :  plastic  state 

(1)  initialize  =  0,  Acr(°l  =  0 

(2)  perform  local  interation 

(3)  update  AA,  Act,  crt+At, 

(4)  check  convergence,  if  not  convergenced  go  to  (2) 

4.  else:  elastic  state,  <Tt+At  =  cr*’" 

Table  3.1;  Updating  the  stress  and  strain  for  Perzyna  model 

3.1.3  Duvaut-Lions  viscoplastic  model 

The  viscoplastic  model  proposed  by  Duvaut  and  Lions  ([4])  is  in  its  elaboration  closer  connected  to  rate 
independent  plasticity.  This  model  is  based  on  the  difference  in  response  between  the  rate  independent 
material  and  the  viscoplastic  material.  Thus  the  model  has  the  advantage  that  it  can  be  combined  with 
a  yield  surface  which  has  an  apex  or  which  is  non-smooth.  The  viscoplastic  strain  rate  and  hardening 


rate  are  defined  as 

€pp=  i[D®]-^(cr-crp) 

T 

(3.28) 

k  =  —  —  {k  —  Kp) 

(3.29) 

where  subscript  p  denotes  the  contribution  of  the  rate  independent  material  and  r  is  the  relaxation  time 
of  the  material  and,  in  general,  is  strain  and  strain  rate  dependent. 

The  stain  rate  decomposition  (i.e.  Eq.  3.1  and  Eq.  3.2)  is  still  used. 

A  two-step  integration  scheme  is  used  to  determine  the  quantities  Cp,  Kp,  a  and  consistent  tangent 
stiffness  matrix  (cf.  [11,  15]). 

-  STEPl;  we  apply  the  fully  implicit  Euler  Backward  method  introduced  in  section  3.1.1  to  integrate 

the  backbone  model  to  determine  Cp  ,  Hp  and  the  tangential  stiffness  relation  (cf.  3.19) 

Aa-p  =  D"PAe  (3.30) 

-  STEP2:  A  generalized  Euler  method  is  used  to  estimate  the  incremental  viscoplastic  strain 

A6i,p  =  [cup(l  —  6)  +  0e‘p'^‘]  At  (3.31) 

where  6  is  the  interpolation  parameter  for  which  0  <  0  <  1. 

Utilizing  the  incremental  form  relations  Eq.  3.10,  Eq.  3.11  and  Eq.  3.28  for  e[,p'^‘  stress  update  and 
material  tangent  stiffness  matrix  can  be  obtained 


Aa  =  D^PAe  -  Aq 


(3.32) 
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1.  update  strain  increment  and  total  strain  Ae,  e 

2.  compute  trial  stress  c*’"  =  D®(e  -  e‘p) 

3.  if  f(cr‘’',K‘)  >  0  :  plastic  state 

(1)  calculate 

(2)  calculate  Atr,  cr*+^\  cr‘+^‘,  4+^‘,  B^p 

4.  if  f(o■‘^  K*)  <  0:  =  cr*'- 


Table  3.2:  Updating  the  stress  and  strain  for  Duvaut-Lions  model 


B^P 


Aq  = 


T  +  9  At 
rAt 

T  +  6  At 


D®  +  —B^P 

T 


{1-9)B%^  +  -<t\ 


(3.33) 

(3.34) 


where  Aq  is  an  extra  pseudo-nodal  force  in  equation  of  motion.  cr(,p  —  is  the  viscous  stress  at 

the  beginning  of  the  time  step. 

In  Table  3.2  the  algorithm  for  Duvaut-Lions  model  is  outlined.  The  update  of  k  at  time  t-fAt  can  be 
performed  by  using 

+  (1  -  e-^*/^)Kp  (3.35) 


Note  that  in  the  Duvaut-Lions  model,  the  loading  and  unloading  condition  is  determined  by  the  total 
stress  state.  In  case  unloading  happens  in  the  back-bone  stress  update  while  the  total  stress  state  is 
still  outside  the  inviscid  yield  surface,  a  modification  should  be  adopted.  The  Kuhn-Tucker  conditions  in 
back-bone  model  (i.e.  Eq.  3.4)  are  modified  as 

/  <  0,  fX  =  0  (3.36) 

A  negative  value  of  A  is  allowed. 

The  maximum  principal  stress  model  of  Rankine  is  used  for  Duvaut-Lions  viscoplasticity 

f{cr,  k)  =  Oi-  &[k)  (3.37) 

in  which  ai  =  max(CTi, (72,0-3).  Considering  a  plane  stress  condition,  we  have 


/(o’:  «)  =  +  Cfy)  4-  \l \{U^  +  Uyf  +  T^y  -  <7(/t) 


(3.38) 


There  is  an  apex  in  the  Rankine  yield  contour,  which  has  been  treated  by  Feenstra,  Pamin  and  Meschke 
(cf.  [5,  8,  7]).  Due  to  the  used  loading  condition  (i.e.  one-directional  monotonic  loading)  the  apex  issue 
is  not  treated  in  the  simulation. 


3.2  Continuum  damage  model 

The  continuum  damage  model  describes  changes  at  the  microstructual  level  in  the  material  via  a  finite 
number  of  scalar  or  tensor-valued  internal  variables.  In  this  report,  an  isotropic  elasticity  based  damage 
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model  with  a  single  scalar-valued  damage  variable  is  used, 

<r=(l-w)D"e  ,  0<a;<l  (3.39) 

where  w  is  a  scalar  damage  parameter.  It  is  a  function  of  history  parameters  and  its  evolution  is  controlled 
by  a  loading  function. 

In  the  local  formulation  of  the  damage  model,  the  damage  parameter  w,  history  parameter  k  and 
loading  function  can  be  defined  as 

u  =  w{k)  (3.40) 

K  =  max{eeq,  K-i)  (3.41) 

/(Ce,,  k)  =  Ceq  -  K  (3-42) 

where  the  equivalent  strain  egq  is  an  invariant  measure  of  strain,  /Cj  is  the  threshold  of  damage  initiation. 

Similar  to  the  standard  plasticity  theory  the  loading-unloading  condition  can  be  formulated  as  Kuhn- 
Tucker  conditions 

/<0,  K>0,  fk  =  0  (3.43) 

Just  as  the  standard  plastic  model  combined  with  strain  softening  law,  at  a  certain  stage  of  damage 
evolution,  the  continuum  damage  model  ceases  to  produce  meaningful  results.  The  governing  equations 
of  motion  become  elliptic  and  the  IVP  becomes  ill-posed. 

Nonlocal  and  gradient  enhancedment  have  the  regularizing  effect  to  the  continuum  damage  model.  The 
constant  gradient  damage  (CGD)  model  and  the  strain-based  transient-gradient  damage  model  (STGD) 
will  be  adopted  in  this  report. 


3.2.1  Implicit  constant  gradient  damage  model 

In  the  nonlocal  formulation  of  the  constitutive  model  the  nonlocal  conterpart  x  of  a  local  state  variable 
X  over  the  surrounding  domain  G  with  respect  to  the  weight  function  G  is  given  by  (cf.  [9,  2,  1]) 


^G(0x(x  +  g)df2(0 


(3.44) 


with  $  a  vector  pointing  to  the  infinitesimal  volume. 

Nonlocal  constitutive  relations  can  be  considered  as  a  point  of  departure  for  constructing  a  gradient 
model.  The  local  quantity  x(x)  can  be  developed  into  a  Taylor  series  around  x  and  substituted  into 
Eq.  3.44,  under  the  assumption  of  isotropy,  a  new  expresssion  for  x  is  found 


X(x)  =  x(x)  +  C'2V2x(x)  +  C4V'x(x)  +  . . . 


(3.45) 


where  V  =  [d/dx,  d/dy,  djdzY ■  This  series  is  called  the  explicit  gradient  enhancement. 

An  alternative  gradient  formulation  can  be  derived  from  Eq.  3.45  by  applying  the  Laplacian  operator 
to  it  and  multiplying  by  C^-  The  result  is  substracted  from  Eq.  3.45  ,  which  gives 

X(x)  -  C'2V2x(x)  =  x(x)  +  (C4  -  C'22)VV(x)  +  (Ce  -  C2C4)V«x(x)  +  . .  •  (3.46) 


This  series  is  called  implicit-gradient  enhancement.  The  same  procedure  as  above  can  be  repeated  to 
derive  higher-order  implicit  gradient-enhancement. 

The  parameters  . . .  depend  on  the  weight  function  G. 

It  has  been  proved  that  implicit  formulations  show  better  performance  than  explicit  formulations. 
Considering  the  equivalent  strain  eeg(x)  as  the  local  state  variable,  the  second-order  implicit  gradient- 
enhancement  can  be  formulated  as 


£eq(x)  -  cV^£eg(x)  =  eeq(x) 


(3.47) 


where  the  gradient  parameter  c  has  the  dimension  of  length  squared.  We  adopt  the  phenomenological 
view  that  c  reflects  the  length  scale  of  the  failure  process  that  we  wish  to  describe  macroscopically. 

Higher-order  continuum  requires  additional  boundary  conditions,  either  of  the  Dirichlet  {leg  —  e)  or  of 
the  Neumann  type  (n^Vce,  =  e„).  The  natural  boundary  condition 

n^Vce,  =  0  (3.48) 

where  n  is  the  outward  normal  to  the  boundary,  is  used  in  the  simulations. 


Now  the  gradient  enhanced  damage  constitutive  relation  can  be  written  as 
tr=(l-a;)D"e  ,  0  <  a;  <  1 


(3.49) 


oj  =  u){k)  (3.50) 

K  =  max{eeq,  Ki)  (3.51) 

/(Ceg,  k)  =  eeg-  k  (3.52) 

/<0,  «>0,  fk=0  (3.53) 


The  nonlocal  equivalent  strain  ieq  is  defined  through  Eq.  3.47 
We  adopt  two  damage  evolution  laws,  namely 

-  the  linear  softening  damage  evolution  law 

=  1  -  (3,54) 

K  Kc  —  Ki 

characterized  in  the  one-dimensional  situation  by  a  linear  decrease  of  stress  until  a  zero-stress  level  is 
reached  for  the  ultimate  strain  Kc, 

-  the  exponential  softening  damage  evolution  law 

a;(/c)  =  1  _  ^  J(1  _  a)  +  (3.55) 

where  a  and  /3  are  constant  parameters, 


and  the  modified  Von  Mises  definition  of  equivalent  strain 

{k-l)h  1  /(fc-l)2/i2  2k  J2 
2k{l-2u)  2k\l  (l-2:/)2  (H- !/)2 


(3.56) 


where  Ii  the  first  invariant  of  strain  tensor,  J2  the  second  invariant  of  the  deviatoric  strain  tensor  and  k 
the  ratio  between  uniaxial  compression  yield  stress  and  uniaxial  tension  yield  stress. 

In  section  2.2  we  have  derived  a  weak  form  of  the  initial  value  problem.  The  (I-l-l)-th  iteration  form 
at  time  t-l-At  is 


f  ]dn -1-  f  6€^dcr^^^dQ=  [  6u^p[!^^^dQ, 

Jn  Jn  Jn 

+  f  -  f  Se^a[^_^^^dn 

JdQt  Jn 


(3.57) 


Similarly,  the  weak  form  of  the  modified  Helmholtz  equation  Eq.  3.47  can  be  formed 
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Using  Green’s  formula  of  integration  by  parts  and  considering  boundary  condition  Eq.  3.48,  the  weak 
form  of  the  diffusion  equation  results  in 


[  [  V5eJ,cVeegdfi  =  f 

Jn  Jci  Jo. 

The  (I+l)-th  iteration  form  at  time  t+At  can  be  written  as 


(3.59) 


S^egei^glt+Atd^ 


-  +  V56%cVe-W  dO 


(3.60) 


where  we  have  introduced  the  decompositions 


?(^+i)  =  g(U 


eg, t+At  ^cg,t+Af  ~  ^^€g,t+At 


+  de 


(U 


^69, t+At  ^eg.t+At  ^  t+At 

Quantities  dcrU)  and  d€eq^^'>  can  be  expressed  as 


.(U 


(U 


da^^^  =  (1  - 


5u 


(U 


dn 

d€eq 


in 


= 


f  deeq^^ 

T 

'de' 

[de  ) 

du(^) 


(3.61) 

(3.62) 

(3.63) 

(3.64) 


where  subscript  t+At  is  dropped  for  simplicity.  In  Eq.  3.63  ,  =  1  for  loading  and  =  0  otherwise. 


For  the  finite  element  implementation,  a  two-field  discretization  is  used 
u(x,t)  =  H„(x)a(t) 

u(x,t)  =  H,,(x)a(t) 

e(x,t)  =  B„(x)a(t)  =  LHua(t) 


(3.65) 

(3.66) 

(3.67) 


£e,(x,t)-He(x)A(t)  (3-68) 

Veeg(x,t)  =  BeA(t)  =  VH,(x)A(t)  (3.69) 

where  a(t)  and  A(t)  contain  the  assembled  nodal  values  of  the  nodal  unknowns,  the  displacements  and 
the  nonlocal  equivalent  strains. 

Substitution  of  Eqs.  3.63  to  3.69  into  the  weak  forms  of  the  field  equations  3.57  and  3.60  yields 
da^M(^)a(^+i)  +  da^K^da^^)  +  da^K^^.^dA^^)  =  da^fi^+P  -  da^C,.  (3-70) 

dA^K(l)da(^)  +  dA7Klf)dA(')  =  -dA^f^^  ^  (3.71) 

where 

M(^)  =  /  HjRHudn  (3.72) 

Jo 

f  B7(l-c^U))D«B„dG 
Jo 


(3.73) 
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K«=  [  (HfH,  +  BjcB,)dn 
Jn 


+  [ 

Jn  Jan 

€lu=  f  BW^Un 

Jn 


Budn 


(3.74) 

(3.75) 

(3.76) 

(3.77) 

(3.78) 


d  =  ^  (Hf  +  B^cB,A(^)  -  Hfe^)  dQ  (3.79) 

Since  equations  3.70  and  3.71  must  hold  for  any  arbitrary  <5a  and  dA  the  semi-discretized  weak  forms 
of  field  equations  then  give 

+  KWda(^)  +  "  €lu  (3-80) 


KWda(')  -f  (3.81) 

We  apply  the  implicit  time  integration  scheme  of  Newmark  family  to  Eqs.  3.80  and  3.81  to  obtain 
fully  discretized  system  of  equations  for  the  initial  value  problem.  Following  the  procedure  described  in 
section  2.4,  we  obtain 


J\.U£ 

da^^^ 

f(/+i)  _  j.(/)  _ 

_f(/) 

(3.82) 


which  can  be  solved  by  using  a  full  Newton-Raphson  procedure.  The  algorithmic  treatment  for  the  global 
iteration  process  is  similar  as  the  one  outlined  in  Table  2.1  in  section  2.4.  While  the  procedure  to  update 
the  stress  and  strain  relation  is  outlined  in  Table  3.3. 


3.2.2  Strain-based  transient-gradient  damage  model 

Enhanced  continuum  models  perform  worse  at  the  later  stages  of  failure  because  they  can  not  reflect 
the  development  of  discrete  surfaces.  Where  a  discrete  surface  should  develop  the  strain  and  strain  rate 
approach  infinity,  which  can  lead  to  spurious  spreading  of  inelastic  deformations,  damage  growth  and  for 
rate  dependent  models  even  a  hardening  effect  since  strain  rate  becomes  very  high  (cf.  [6,  16]). 

In  order  to  represent  the  complete  failure  process,  displacement  discontinuities  can  be  incorporated  at 
certain  stages  or  remeshing  techniques  should  be  used.  For  the  implicit  gradient  damage  model  described 
in  the  previous  section,  Geers  ([6])  proposed  to  adopt  a  variable  length  scale,  and  eliminate  the  interaction 
between  crack  and  the  surrounding  material. 

The  governing  system  of  partial  differential  equations  now  read  (cf.  Eq.  2.1  and  Eq.  3.47) 

L^cr  -1-  p  =  Rii  in  0,  (3.83) 

£e5(x)  -  CV^£eg(x)  =  ee,(x)  in  Q  (3.84) 

Instead  of  a  constant  gradient  parameter  as  in  Eq.  3.47,  a  variable  C  is  used  which  is  denoted  as  gradient 
activity.  It  models  the  mobilized  nonlocal  coupling  between  particles  at  microlevel. 
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1.  update  strain  increment,  total  strain  and  nonlocal  equivalent  strain  at  integration  point 

Ac  ,  e  ,  €eg 

2.  compute  equivalent  strain  eeg(^)> 

3.  evaluate  the  loading  function  f  = 

where  kq  is  the  converged  value  of  k  at  the  end  of  the  previous  load  increment 

4.  if  loading(f  >  0)  then  k  =  Cg, 

else  K  =  Ko 

5.  compute  damage  variable  and  its  derivative  with  respect  to  cu  =  u>(k), 

6.  compute  the  new  stress  tr  =  (1  -  u;)D®e 

Table  3.3:  Updating  the  stress  and  strain  for  implicit  gradient  damage  model 


de. 

^eg  "  ^0 


For  the  strain-based  transient-gradient  damage  model,  we  use  the  following  gradient  activity  evolution 
law 

^^^eq<£eqX  (3.85) 

I.  C  if  Ceq  t>  ^eqX 

where  c  and  are  two  constants.  This  evolution  law  implies  that  nonlocal  coupling  is  mobilized  more  and 
more  as  the  local  deformation  increases,  while  it  gradually  vanishes  in  the  unloaded  material  surrounding 
cracks.  Thus  the  large  strain  in  cracks  gives  little  effect  to  surrounding  material,  and  the  spurious  damage 
growth  at  the  crack  faces  is  prevented. 

A  three-field  formulation  is  used  to  achieve  stability,  [u,ee^,c].  The  weak  forms  of  the  governing 
equations  are  written  as  (cf.  Eq.  2.11  and  Eq.  3.58) 


f  (5u^[L^ct  +  p  —  Rujdfi  =  0 

Jn 

Jn  Jn 


(3.86) 

(3.87) 


/ 

Jn 


5<r(<  -  c)dn  =  0 


(3.88) 


For  implementation  using  the  finite  element  method,  the  continuous  fields  u(x,t),  eeg(x,t)  and  <r(x,t)  are 
discretized  by  using  shape  functions. 


u(x,  t)  =  Hu(x)a(t) 


(3.89) 


u(x,t)  =  H„(x)a(t) 


(3.90) 


e(x,t)  =  Bu(x)a(t)  =  LHua(t) 


(3.91) 
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eeq{x,t)  =  H£(x)A(i) 


V6e,(x,  t)  =  B,A(t)  =  VH,(x)A(t)  (3.93) 

?(x,i)  =  H^(x)T(f)  (3.94) 

V<r(x,i)  =  B,T(i)  -  VH,(x)T(i)  (3.95) 

With  the  aid  of  Green’s  formula  of  integration  by  parts,  we  obtain 

f  5u^[Ru]dn  4-  [  Se^udQ  =  f  du^pdQ  +  f  Su^tdT  ‘  (3.96) 

Ja  Jn  Jn  Jaut 

[  [Se^gSeq  +  VSiJ  <;Veeq  +  5€eqV<;'^Veeq]  dQ=  I  SeeqesqdCl  (3.97) 

Jn  Jn 

[  (5<r[?-C]dn  =  0  (3.98) 

Jn 

where  boundary  condition  N’^cr  =  t  (i.e.  Eq.  2.5)  and  additional  boundary  condition  qVeeq  n  —  0  (cf. 
Eq.  3.48)  are  used. 

A  linearization  process  and  time  integration  scheme  similar  to  that  described  in  section  3.2.1  are  used. 
Then  a  full  discrete  system  of  governing  equations  is  obtained  (cf.  3.82) 

■  Ki^J  +  CoMG)  Ki{^  0  1  r  daG)  ]  f  fiitu  -  ' 

kIP  KiP  dAG)  =  -fG)  (3.99) 


0 


where  subscripts  t+At  are  omitted  for  simplicity.  The  expressions  of  the  sub-matrices  and  the  vectors 
are  listed  below 


M('f)  =  f  HjRH^dfl 

Jn 

KV>  = 

kS  =  -  /  BlD-em 


KS  =  -JnT[{^)  \  B,dn 

KW  =  j  HfH,-t-BfH,T(^^B,-pHj  [y(^)]^B^B,  dfl 
=  J  Bj B, +  Hf  ^  B^ B, j  dG 


=  -  [  HjU,dQ 
Jn 

Jn  Jant 


(3.100) 

(3.101) 


(3.102) 


(3.103) 

(3.104) 

(3.105) 

(3.106) 

(3.107) 

(3.108) 

(3.109) 


C.u  =  [  (3.109) 

4m  =  J  +  BfH,T('^)B,A(^)  +  Hf  [T^^^j^BfB.A^-^^dn  -  J  (3.110) 

41,  =  -  /  H^4)dG+  [  H^H,TWdG  (3.111) 

’  Jn  Jn 

The  procedure  for  updating  stress  and  strain  at  integration  point  level  is  the  same  as  in  Table  3.3, 
while  the  global  incremental-iterative  procedure  is  similar  to  that  outlined  in  Table  2.1. 
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Chapter  4 


One-dimensional  localization  problem 


In  this  chapter,  we  treat  a  one-dimensional  bar  problem  under  impact  tensile  loading.  The  continuum 
models  introduced  in  Chapter  3  are  used.  The  main  purpose  is  to  demonstrate  the  numerical  consequence 
of  the  ill-posedness  of  the  initial  value  problem  for  standard  plastic  model  and  local  damage  model,  and 
the  regularization  effect  of  rate  dependent  viscoplasticity  models  and  gradient  enchanced  damage  models. 
Only  the  mode-I  localization  problem  is  considered.  The  results  can  be  extended  to  the  mode-II  problem 
if  an  appropriate  yield  function  is  used.  The  one-dimensional  bar  used  in  the  numerical  simulation  is 
sketched  in  Fig.  4.1.  The  consistent  mass  matrix  is  adopted  in  the  simulation. 


Figure  4.1:  The  one-dimensional  bar  problem. 


4.1  Standard  elastoplastic  model 

The  standard  plastic  model  described  in  section  3.1.1  is  tested  first.  A  rankine  yield  criterion,  an  as¬ 
sociative  flow  rule,  a  linear  softening  law  and  a  strain  hardening/softening  hypothesis  are  utilized.  The 
material  parameters  are  listed  in  Table  4.1.  Use  has  been  made  of  a  two-noded  line  element.  The  bar  is 
divided  into  10,  20,  40  and  80  elements,  respectively.  A  constant  tensile  impact  velocity  V  =  450  mm/s  is 
applied  on  the  loading  end  of  the  bar,  the  bar  responses  linear  elastically  until  the  loading  wave  reaches 
the  left  boundary.  The  doubling  in  stress  due  to  reflection  of  tensile  wave  causes  the  initiation  of  yielding. 
The  material  enters  the  softening  regime  and  strain  localization  zone  emerges.  A  time  step  At  =  0.5  fj,s 
is  used  for  all  calculations. 

Table  4.1:  Material  parameters  for  standard  elastoplastic  model  in  the  one-dimensional  problem. 


E 

En 

cro 

P 

3.5e-f4  MPa 

-3.5e-|-3  MPa 

5.0  MPa 

2.0e-9  Ns^/mnP 

The  strain  and  stress  profiles  for  different  meshes  at  time  t  =  36.0  /is  w  (c®  is  elastic  wave 
velocity)  are  shown  in  Fig.  4.2.  It  can  be  seen  that  strain  localization  occurs  in  one  element  (one 
integration  point)  which  is  the  smallest  possible  zone.  Hence  the  width  of  the  localization  zone  decreases 
upon  mesh  refinement.  The  amount  of  wave  reflection  depends  on  the  discretization.  As  soon  as  the 
stress  has  become  zero,  the  localization  zone  (i.e.  one  integration  point  near  the  left  boundary)  acts  as 
a  free  boundary  on  which  tensile  wave  reflects  as  a  pressure  wave  (cf.  [11]).  Note  that  all  results  for  the 
one-dimensional  problem  are  the  values  of  variables  at  individual  integration  points. 
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20 


(a)  Strain  (b)  Stress 

Figure  4.2:  The  strain  profiles  (a)  and  stress  profiles  (b)  for  different  meshes. 


4.2  Perzyna  viscoplastic  model 

The  Rankine  yield  criterion  combined  with  Perzyna  viscoplasticity,  an  associative  flow  rule,  a  linear 
softening  law  and  a  strain  hardening/softening  hypothesis  are  used.  The  material  parameters  are  given 
in  Table  4.2.  Use  has  been  made  of  a  two-noded  line  element.  A  constant  impact  tensile  velocity  V  = 
450  mm/s  is  applied  on  the  loading  end. 


Table  4.2:  Material  parameters  for  Perzyna  viscoplastic  model  in  one-dimensional  problem . 


E 

En 

(^0 

V 

N 

3.5e+4  MPa 

-3.5e+3  MPa 

5.0  MPa 

1.0 

(b)  Stress 


(c)  Internal  length  scale 


Figure  4.3:  The  strain  profiles  for  different  meshes  (a),  the  stress  profiles  for  different  meshes  (b)  and  the 
strain  profiles  for  different  internal  length  scales  (c). 


The  strain  and  stress  profiles  for  different  discretizations  (20,  40,  80,  160  and  320  elements)  at  time 
t  =  36.0  /iS  are  shown  in  Figs.  4.3a  and  b.  Constant  time  steps  At  =  0.10  ^is  for  the  calculation  with 
320  elements  and  At  =  0.15  fj.s  for  other  calculations  are  used.  The  results  show  an  exponential  decrease 
of  strain  after  the  reflection  on  the  fixed  end.  For  the  coarse  meshes  of  20,  40  and  80  elements,  the 
results  are  not  identical.  However,  a  unique  solution  is  obtained  when  the  mesh  is  refined  to  160  and 
320  elements.  The  localization  zone  converges  to  a  finite,  constant  band  width  upon  mesh  refinement. 
The  mesh  independent  partial  reflection  on  the  localization  zone  can  be  seen  from  the  stress  profiles  in 
Fig.  4.3b.  This  result  shows  that  the  viscous  effect  has  a  regularization  effect. 

The  influence  of  viscosity  parameter  rj  on  the  width  of  localization  zone  can  be  seen  in  Fig.  4.3c  for 
the  mesh  with  160  elements.  The  width  of  localization  zone  increases  with  the  decrease  of  viscosity 
parameter,  i.e.  the  increase  of  viscous  effect. 
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4.3  Duvaut-Lions  viscoplastic  model 

The  Rankine  yield  criterion  combined  with  Duvaut-Lions  viscoplasticity,  an  associative  flow  rule,  a  linear 
softening  law  and  a  strain  hardening/softening  hypothesis  are  used.  The  material  parameters  are  given 
in  Table  4.3.  Use  has  been  made  of  a  two-noded  line  element.  A  constant  impact  tensile  velocity  V  = 
450  mm/s  is  applied  on  the  loading  end. 


Table  4.3:  Material  parameters  for  Duvaut-Lions  viscoplastic  model  in  one-dimensional  problem. 


E 

En 

ffo 

r 

P 

3.5e-l-4  MPa 

-3.5e-l-3  MPa 

5.0  MPa 

1.0  /xs 

2.0e-9  Ns^/mrrP 

(c)  Internal  length  scale 


Figure  4.4:  The  strain  profiles  for  different  meshes  (a),  the  stress  profiles  for  different  meshes  (b)  and  the 
strain  profiles  for  different  internal  length  scales  (c). 

The  strain  and  stress  profiles  for  different  discretizations  (20,  40,  80,  160  and  320  elements)  at  time  t 
=  36.0  fj.s  are  shown  in  Figs.  4.4a  and  b.  Constant  time  steps  At  =  0.10  /xs  for  the  calculation  with  320 
elements  and  At  =  0.15  /xs  for  other  calculations  are  used.  The  convergence  to  a  unique  solution  upon 
mesh  refinement  for  the  strain  localization  and  wave  reflection  show  that  physically  realistic  results  can 
be  obtained.  The  influence  of  the  relaxation  time  r  on  the  width  of  the  localization  zone  is  presented  in 
Fig.  4.4c. 

For  the  one-dimensional  case,  a  dispersion  analysis  can  be  carried  out  for  the  Perzyna  and  the  Duvaut- 
Lions  viscoplastic  model  (cf.  [11]).  The  obtained  internal  length  scale  for  the  Perzyna  model  is  1  = 
for  mode-I  localization.  A  similar  expression  for  the  Duvaut-Lions  model  is  1  =  2rCe-  When  the  internal 
length  scales  1  for  the  two  models  are  equal,  almost  identical  numerical  results  can  be  obtained. 


4.4  Local  damage  model 

A  local  damage  model  combined  with  a  linear  softening  damage  evolution  law  and  an  equivalent  strain 
eeq  =  €  is  adopted  in  this  section.  The  material  parameters  are  presented  in  Table  4.4  .  Use  has  been 
made  of  a  two-noded  line  element.  The  bar  is  divided  into  10,  20,  40,  80  elements.  A  constant  impact 
tensile  velocity  V  =  450  mm/s  is  applied  on  the  loading  end.  A  time  step  At  =  0.5  /xs  is  used  for  all 
calculations. 


Table  4.4:  Material  parameters  for  local  damage  model  in  one-dimensional  problem. 


E 

Ki 

Kc 

P 

3.5e-t-4  MPa 

1.43e-4 

1.43e-3 

2.0e-9  Ns^/mrrP 

The  strain,  stress  and  damage  profiles  at  time  t  =  36.0  /xs  for  different  meshes  are  shown  in  Fig.  4.5. 
Strain  localization  occurs  in  one  element  which  is  the  smallest  possible  zone.  The  width  of  the  strain 
localization  zone  and  the  amplitude  of  strain  and  damage  exhibit  mesh  dependence  (cf.  [2]  ).  The  results 
are  similar  to  that  of  a  standard  elastoplastic  model  in  section  4.1. 
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(c)  Damage 


Figure  4.5:  The  strain  profiles  (a),  stress  profiles  (b)  and  damage  profiles  (c)  for  different  meshes. 


4.5  Constant  gradient  damage  model 

The  second-order  implicit  constant  gradient  damage  model,  described  in  section  3.2.1,  is  used  to  illustrate 
the  reguralizing  effect  of  the  gradient  enhancement.  The  linear  softening  damage  evolution  law  and  the 
equivalent  strain  €eq  =  e  are  applied  in  this  model.  Use  has  been  made  of  a  three-noded  line  element  with 
quadratic  interpolation  for  the  displacement  and  linear  interpolation  for  the  nonlocal  equivalent  strain 
ieq.  A  two-point  Gauss  quadrature  scheme  is  used.  The  material  parameters  are  listed  in  Table  4.5.  A 
constant  impact  tensile  velocity  V  =  450  mm/s  is  applied  on  the  loading  end. 


Table  4.5:  Material  parameters  for  constant  gradient  damage  model  in  one-dimensional  problem. 


E 

Ki 

3.5e-l-4  MPa 

1.43e-4 

1.43e-3 

Figure  4.6:  The  strain  profiles  (a),  stress  profiles  (b)  and  damage  profiles  (c)  for  different  meshes. 

The  numerical  results  are  presented  in  Fig.  4.6  including  the  strain,  stress  and  damage  profiles  for 
different  discretizations  (20,  40,  80,  160  and  320  elements)  at  time  t  =  36.0  /us.  Constant  time  steps 
At  =  0.10  ^s  for  the  calculation  with  320  elements  and  At  =  0.15  fis  for  other  calculations  are  used. 
Upon  mesh  refinement,  the  numerical  solution  shows  convergence  towards  a  unique  solution.  Whereas 
the  20  element  mesh  is  still  too  coarse,  the  40  element  mesh  already  gives  a  reasonable  prediction  of 
the  width  of  the  localization  zone.  Hence  the  implicit  gradient  enhanced  damage  model  can  carry  out 
failure  analysis  in  an  objective  manner.  However,  the  mesh  must  be  fine  enough  in  order  to  capture  the 
intense  straining  in  the  localization  zone.  This  holds  true  for  a  higher-order  spatial  derivatives  enhanced 
continuum  model  as  well  as  for  a  higher-order  temporal  derivatives  enhanced  continuum  model.  The 
influence  of  gradient  parameter  c  on  the  width  of  localization  zone  is  presented  in  Fig.  4.7.  The  width 
of  localization  zone  increases  with  the  increase  of  the  parameter  c.  In  the  calculations  shown  in  this 
figure,  the  mesh  of  160  elements  is  used.  The  stress  oscillation  is  due  to  artificial  higher-order  frequencies 
introduced  by  semi-discretization  of  the  equation  of  motion.  Algorithmic  damping  can  be  used  to  exclude 
these  high-frequent  responses.  However,  by  introducing  numerical  damping  in  the  Newmark  scheme  the 
accuracy  deteriorates. 
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(a)  Strain 


(b)  Damage 


Figure  4.7:  The  strain  profiles  (a)  and  damage  profiles  (b)  for  different  gradient  parameter  c. 


4.6  Strain-based  transient-gradient  damage  model 

Because  of  the  deficiency  of  the  constant  gradient  enhanced  damage  model  at  final  stages  of  failure  a 
variable  gradient  enchanced  damage  model  has  been  introduced  in  section  3.2.2.  In  this  section,  the 
one-dimensional  dynamic  tensile  bar  problem  is  used  to  compare  these  two  models  numerically.  Instead 
of  the  linear  softening  law,  an  exponential  softening  law  is  used.  The  equivalent  strain  is  defined  as  te, 
=  e.  The  material  parameters  for  both  models  are  listed  in  Table  4.6  and  4.7,  respectively.  Use  has 
been  made  of  a  three-noded  line  element,  where  quadratic  interpolation  for  the  displacement  field,  linear 
interpolation  for  the  nonlocal  equivalent  strain  field  and  linear  interpolation  for  the  gradient  activity  field 
c  are  utilized.  The  bar  is  divided  into  160  elements  with  At  =  0.15  /iS  for  the  constant  gradient  damage 
model  and  640  elements  with  At  =  0.0375  ps  for  a  strain-based  transient-gradient  damage  model.  A 
constant  impact  tensile  velocity  V  =  590  mm/s  is  applied  on  the  loading  end  of  the  bar. 


Table  4.6:  Material  parameters  for  constant  gradient  damage  model. 


E 

lii 

a 

P 

3.5e-l-4  MPa 

1.43e-4 

0.96 

2400.0 

Table  4.7:  Material  parameters  for  strain-based  transient-gradient  damage  model. 


E 

Ki 

a 

P 

C 

1 

3.5e-b4  MPa 

1.43e-4 

0.96 

2400.0 

9.0  mmP 

1  1.43e-3 

Figure  4.8:  Stroboscopic  development  of  strain  (a),  stress(b)  and  damage  (c)  for  constant  gradient  damage 
model. 

The  numerical  result  for  a  constant  gradient  damage  model  predicts  a  physically  unacceptable  growth 
of  the  damage  zone  when  the  damage  in  the  localization  zone  approaches  the  critical  value  (i.e.  a;  1). 
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X  (iwn)  *  (nwn)  *  <"«") 

(a)  Strain  (b)  Damage  (c)  Gradient  activity 

Figure  4.9:  Stroboscopic  development  of  strain  (a),  damage  (d)  and  gradient  activity  (c). 

This  is  presented  in  Fig.  4.8.  The  reason  is  that  the  intense  increase  of  local  strain  in  the  localization  zone 
results  in  an  increase  of  the  local  equivalent  strain.  The  diffusion  equation  maps  this  increase  onto  the 
nonlocal  equivalent  strain,  which  provokes  an  increase  of  damage.  The  weakening  of  the  damaged  zone 
results  again  in  larger  strains.  The  calculation  is  conducted  until  t  =  42.0  fis  and  the  interval  between 
the  plots  is  3.0  /xs.  The  damage  plateau  of  the  small  amplitude  in  the  other  part  of  the  bar  is  due  to  the 
stress  oscillation  of  the  loading  wave.  It  can  be  avoided  when  the  time  step  is  reduced  to  At  =  0.10  fis. 

The  stroboscopic  development  of  strain  and  damage  for  the  strain-based  transient-gradient  damage 
model  are  given  in  Figs.  4.9a  and  b.  It  can  be  seen  that  this  model  can  preserve  the  width  of  localization 
zone  and  avoid  the  spurious  spreading  of  damage.  The  development  of  gradient  activity  is  presented  in 
Fig.  4.9c.  Due  to  the  steep  change  in  gradient  activity,  the  mesh  must  be  finer  than  that  used  for  a 
constant  gradient  damage  model. 


Chapter  5 

Two-dimensional  microstructure  analysis 


In  this  chapter,  a  few  square  block  microstructures  with  a  width  of  13.774mm  are  used  as  the  represen¬ 
tative  volume  element  for  concrete  to  investigate  the  dynamic  response  of  concrete  at  meso-level.  The 
continuum  models  described  in  chapter  3  are  used.  The  structures  and  the  finite  element  discretization 
are  shown  in  Fig.  5.1.  Three  microstructures  are  used  with  two  meshes  for  stucture  one.  The  microstruc¬ 
tures  are  composed  of  aggregates,  matrix  and  an  interfacial  transition  zone  (ITZ).  Hence,  three  sets  of 
material  parameters  are  used.  A  three-noded  triangular  element  is  used  for  the  discretization  and  a 
lumped  mass  matrix  is  adopted  for  the  numerical  simulation.  The  plane  stress  condition  is  assumed.  A 
constant  time  step  At  =  0.01  fj,s  is  chosen  for  all  the  calculations.  The  boundary  and  loading  conditions 
are  illustrated  in  Fig.  5.2,  the  bottom  of  the  microstructure  is  fixed  while  at  the  top  edge  a  constant 
impact  tensile  velocity  is  applied. 


(a)  Structure  one:  mesh  lA;  4487  nodes, 
8788  elements 


(c)  Structure  two:  5815  nodes,  11423  ele¬ 
ments 


(b)  Structure  one:  mesh  IB;  7317  nodes, 
14327  elements 


(d)  Structure  three:  7951  nodes,  15640  ele¬ 
ments 


Figure  5.1:  Microstructure  and  finite  element  discretization. 
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Figure  5.2:  Boundary  and  loading  conditions  for  the  microstructures. 


5.1  Standard  elastoplastic  model 

The  material  parameters  are  shown  in  Table  5.1  for  the  standard  elastoplastic  model  described  in  section 
3.1.1.  A  Von  Mises  yield  criterion,  an  associative  flow  rule,  an  exponential  softening  law  and  a  strain 
hardening/softening  hypothesis  are  utilized  in  the  material  model. 


Table  5.1:  Material  parameters  for  standard  elastoplastic  model 


Matrix 

Aggregate 

ITZ 

E  3.5e+4  MPa 

V  0.3 

do  5.25  MPa 

doo  5.25e-2  MPa 

(3  200.0 

p  2.0e-9  Ns'^/mm3 

E  5.5e+4  MPa 

u  0.3 

do  41.25  MPa 

doo  41.25e-2  MPa 

P  200.0 

p  2.0e-9  Ns^ /mmP 

E  3.5e+4  MPa 

1/  0.3 

do  3.5  MPa 
doo  3.5e-2  MPa 

P  200.0 

p  2.0e-9  Ns'^/mmP 

A  constant  impact  tensile  velocity  V  =  300  mm/s  is  applied  at  the  loading  edge  of  microstructure  one. 
The  development  of  the  plastic  strain  measure  k  for  mesh  lA  (cf.  Fig.  5.1a)  is  shown  in  Fig.  5.3.  We 
observe  that  before  the  loading  waves  arrive  at  the  fixed  edge  (see  Fig.  5.3a),  the  material  in  some  ITZs 
has  already  failed.  Plastic  straining  is  localized  in  the  ITZs.  After  the  loading  waves  are  reflected  at  the 
fixed  bottom  edge  the  localization  zones  in  the  ITZs  begin  to  propagate  into  the  matrix.  At  time  t  =  8.0 
/iS,  the  waves  which  were  reflected  from  the  fixed  edge  have  reflected  on  the  loading  edge.  Due  to  this 
further  loading,  the  localization  zones,  initiated  from  the  ITZs,  propagate  through  the  matrix  and  connect 
with  each  other,  while  the  plastic  strain  level  is  comparable  with  that  in  the  ITZs.  No  localization  zones 
can  be  seen  initiating  from  the  matrix  material.  The  aggregates  can  not  fail  during  the  this  process. 


Figure  5.3:  The  development  of  history  parameter  k. 


The  width  of  the  localization  zone  is  about  the  size  of  one  elment.  This  is  exhibited  by  the  deformation 
pattern  plot  in  Fig.  5.4a,  which  is  amplified  with  a  factor  of  1200.  When  a  finer  discretization  mesh  IB 
is  used,  the  localization  zones  are  still  concentrated  in  one  row  of  elements,  see  Fig.  5.4b.  This  is  the 
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(a)  mesh  lA,  t  =  8.0  fis 


(b)  mesh  IB,  t  =  8.0  fjs 


Figure  5.4:  The  deformation  patterns  for  two  discretizations  at  time  t  =  8.0  fis. 


(a)  Central  vertical  cross  section  plot 


(b)  Reaction  stress  history 


Figure  5.5:  The  profiles  of  parameter  k  along  the  central  vertical  cross  section  at  time  t  =  8.0  /as  (a)  and 
the  average  reaction  stress  component  (7yy  versus  time  (b). 


(a)  mesh  lA,  t  =  8.0  /as 


Figure  5.6:  The  deformation  pattern  with  an  impact  velocity  V  =  350  mm/s  (a)  and  the  profiles  for  the 
average  reaction  stress  component  ayy  versus  time  with  different  impact  velocities. 


numerical  consequence  of  the  ill-posedness  of  the  initial  value  problem  when  a  classical  continuum  model 
combined  with  a  softening  material  law  is  utilized.  The  mesh  dependence  of  the  solution  is  also  shown  in 
Fig.  5.5,  where  the  profiles  of  k  along  the  central  vertical  cross  section  of  the  microstructure  at  t  =  8.0 
fjLS  for  different  meshes  are  plotted  (left)  and  the  average  reaction  stress  component  cTyy  along  the  loading 
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edge  versus  time  is  shown  (right). 

Fig.  5.6  presents  the  numerical  results  for  different  impact  velocities.  The  peak  value  of  the  average 
reaction  stress  increases  with  the  increase  of  the  impact  velocity  when  the  impact  velocity  is  relatively 
low.  When  the  impact  velocity  is  further  increased,  a  more  intense  plastic  straining  appears  in  the  matrix 
and  ITZs  and  a  negative  loading  rate  dependence  is  obtained. 


5.2  Perzyna  viscoplastic  model 

For  the  Perzyna  viscoplastic  model,  a  Von  Mises  yield  criterion,  an  associative  plastic  flow  rule,  an  expo¬ 
nential  softening  law  and  a  strain  hardening/softening  hypothesis  are  applied.  The  material  parameters 
are  listed  in  Table  5.2. 


Table  5.2:  Material  parameters  for  Perzyna  viscoplastic  model 


Matrix 

Aggregate 

ITZ 

E  3.5e-t-4  MPa 

V  0.3 

do  5.25  MPa 

doo  0.0 

/?  0.67e-|-3 

p  2.0e-9  Ns^  jmm^ 

N  1.0 

E  5.5e-l-4  MPa 

V  0.3 

do  41.25  MPa 

doo  0.0 

13  1.0e-h3 

p  2.0e-9  Ns^/mrn^ 

N  1.0 

E  3.5e-!-4  MPa 

V  0.3 

do  3.5  MPa 

doo  0.0 

/?  l.Oe+3 

p  2.0e-9  Ns^/mm^ 

N  1.0 

Figure  5.7:  Development  of  plastic  strain  measure  k  for  Perzyna  viscoplastic  model. 

Structure  one  is  analyzed  first.  The  development  of  the  plastic  strain  measure  «;  for  V  =  500  mm/s,  77  = 
200  is  shown  in  Figs.  5.7a  to  c.  Similar  to  that  of  a  standard  plastic  model,  strain  localizations  occur 
in  the  ITZs  before  the  loading  waves  arrive  at  the  fixed  edge.  After  reflection,  two  strain  localization 
zones  initiate  from  the  fixed  corners  of  the  structure.  At  time  t  =  8.0  fis,  the  waves  have  reflected  from 
the  loading  edge,  it  can  be  seen  that  a  few  lower  strain  level  localization  zones  have  initiated  from  the 


0  26-06  46-06  66-06  66-06  0  26-06  46-06  66-06  86-06 

<  (6)  <  (•) 

(a)  Different  impact  velocities  (b)  Different  viscosity  parameters 

Figure  5-9:  Influence  of  impact  velocity  and  viscosity  parameter  on  the  reaction  stress  history  along  the 
loading  edge  for  Perzyna  viscoplastic  model. 

loading  edge  while  the  localization  zones  formed  from  the  fixed  corners  were  constraint  around  those 
corners.  These  localization  zones  in  the  matrix  no  longer  propagate  in  one  row  of  element  and  each  zone 
is  spread  over  a  few  elements  in  its  transverse  direction.  The  plastic  straining  is  mainly  concentrated  in 
the  ITZs,  the  localization  zones  which  are  formed  in  the  ITZs  are  unlikely  to  propagate  into  the  matrix. 
Increasing  the  impact  velocity  to  V  =  600  mm/s  shows  that  there  are  a  few  more  localization  zones  being 
formed  in  the  matrix  after  the  loading  waves  have  reflected  on  the  fixed  edge  and  they  can  connect  with 
each  other,  but  the  most  intense  plastic  straining  is  still  located  in  the  ITZ  (see  Fig.  5.7d  and  Fig.  5.8b). 

This  result  compares  well  with  the  analysis  conducted  by  Weimin  Wang  (cf.  [15]),  which  shows  that 
the  localization  zone  may  be  determined  by  the  imperfection  size  or  the  internal  length  scale  introduced 
by  viscous  effect  depending  on  which  one  is  smaller.  In  the  two-dimensional  case,  the  influence  of  the 
imperfection  is  very  significant  near  the  imperfection  zone,  while  the  internal  length  scale  tends  to  be 
dominant  at  some  distance  away  from  the  imperfection.  Obviously,  because  of  a  lower  strength,  ITZs  act 


(a)  Profiles  of  k  along  the  central  vertical 
cross  section  at  time  t  =  8.0  /rs 


(b)  History  of  average  reaction  stress  com¬ 
ponent  ayy  along  the  loading  edge 


Figure  5.10:  Mesh  objectivity  for  the  Perzyna  viscoplastic  model  with  V  =  600  mm/s,  rj  —  200  s  ^ 


(a)  Deformation  pattern  0)  Stroboscopic  plot 


Figure  5.11:  The  deformation  pattern  at  time  t  =  8.0  fis  for  mesh  lA  (a)  and  stroboscopic  development 
of  plastic  strain  measure  k  along  the  central  vertical  cross  section  for  mesh  IB  (b)  with  V  =  600  mm/s 
and  rj  =  400  ^  I  ^  I  . 


(d)  t  =  2.0  /is  (e)  t  =  5.5  ps  (f)  t  =  8.0  ps 


Figure  5.12:  Development  of  plastic  strain  measure  k  for  structure  two  (a)-(c)  and  structure  three  (d)-(f) 
with  V  =  600  mm/s  and  rj  =  400  s"^. 
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(b)  Structure  three 


Figure  5.13:  Stroboscopic  development  of  plastic  strain  measure  k  along  the  central  vertical  cross  section 
for  structure  two  (a)  and  structure  three  (b);  The  reaction  stress  component  dyy  along  the  loading  edge 
for  the  three  structures  (c). 


(a)  Structure  two  (b)  Structure  three 

Figure  5.14:  The  deformation  pattern  at  time  t  =  8.0  ^s  for  structure  two  (a)  and  three  (b)  with  V  = 
600  mm/s  and  rj  =  400  s“^. 


as  imperfections  in  the  microstructure  of  concrete.  The  aggregates  in  concrete  have  a  high  strength  and 
elastic  modulus.  On  one  hand,  the  partial  reflection  of  loading  waves  on  them  results  in  earlier  and  more 
intense  plastic  straining  in  the  ITZ.  On  the  other  hand,  they  obstruct  the  propagation  of  a  localization 
zone  formed  in  the  matrix.  The  fact  that  some  localization  zones  in  the  matrix  are  concentrated  around 
the  corners  of  the  fixed  edge  (see  Fig.  5.7d)  is  an  example  of  the  latter  role. 

The  influence  of  the  impact  velocity  and  the  viscosity  parameter  is  given  in  Fig.  5.8.  The  results  show 
that  the  plastic  straining  is  more  intense  when  the  impact  velocity  or  viscosity  parameter  are  increased. 
From  the  average  reaction  stress  component  ayy  along  the  loading  edge  versus  time  plots  (see  Fig.  5.9), 
we  observe  that  the  maximum  reaction  stress  level  increases  apparently  with  the  increase  of  the  impact 
velocity  and  the  decrease  of  the  viscosity  parameter  (which  corresponds  to  an  increase  of  the  viscous 
effect.) 

Mesh  objectivity  of  the  numerical  solution  is  checked  in  Fig.  5.10  with  V  =  600  mm/s  and  rj  =  200 
s“^.  It  shows  that  a  finer  mesh  is  needed  in  the  ITZs  in  order  to  capture  the  high  strain  gradient. 

The  deformation  pattern  amplified  with  a  factor  of  1200  is  presented  for  V  =  600  mm/s  and  r/  = 
400  at  time  t  =  8.0  fis  in  Fig.  5.11a,  where  the  discretization  mesh  lA  is  used.  Large  deformation 
mainly  occurs  in  the  ITZs  of  the  microstructure.  The  width  of  the  ITZ  is  very  small,  so  that  the  higher 
strain  level  in  these  regions  does  not  result  in  significant  displacement  change.  Fig.  5.11b  presents  the 
stroboscopic  plot  of  development  of  the  plastic  strain  measure  k  along  the  central  vertical  cross  section 
of  the  microstructure  for  the  same  impact  velocity  and  viscosity  parameter.  As  in  all  stroboscopic  plots 
presented  in  this  chapter,  the  interval  between  the  first  three  profiles  is  1.0  /iS  and  0.5  ^s  between  other 
profiles.  Comparing  with  Fig.  5.1b,  it  can  be  seen  that  strain  localization  in  the  matrix  occurs  mainly 
after  the  loading  waves  have  propagated  back  from  the  fixed  edge  and  have  reflected  on  the  loading  edge. 
It  is  noted  that  the  strain  level  in  the  microstructure  is  not  very  high. 
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Figure  5.15:  Development  of  plastic  strain  measure  k  in  microstructure  one  with  mesh  lA  for  V  =  400 
mm/s  (a)-(c),  V  =  600  mm/s  (d)-(f)  and  t]  =  200  s~F 


(c)  Reaction  stress 


Figure  5.16:  Stroboscopic  development  of  plastic  strain  measure  k  along  the  vertical  line  across  the  center 
of  the  aggregate  at  the  lower  left  corner  of  microstructure  one  with  t?  =  200  The  reaction  stress 
component  ayy  along  the  loading  edge  for  the  two  computations. 


For  the  microstructures  two  and  three  the  impact  tensile  velocity  V  =  600  mm/s  and  the  viscosity 
parameter  rj  =  400  s“^  are  utilized.  The  development  of  plastic  strain  measure  k  is  shown  in  Fig.  5.12.  As 
in  structure  one  the  strain  localization  is  concentrated  in  the  ITZs.  The  localization  zones  initiated  from 
the  fixed  corners  of  the  structure  are  confined  by  the  aggregates  and  can  not  develop  into  a  dominant  shear 
band.  In  microstructure  three  a  few  more  aggregates  are  placed  around  one  of  the  two  fixed  corners,  it 
can  be  seen  that  the  localization  zone  can  not  initiate  in  the  matrix  around  that  corner.  The  stroboscopic 
development  of  k  along  the  central  vertical  cross  section  for  these  examples  are  presented  in  Figs.  5.13a 
and  b.  The  deformation  patterns  are  given  in  Fig.  5.14. 

The  plot  in  Fig.  5.13c  compares  the  average  reaction  stress  component  eXyy  along  the  loading  edge 
for  the  three  structures  with  mesh  lA  for  structure  one.  We  can  observe  that  the  peak  stress  level 
decreases  with  the  increase  of  the  volume  of  aggregates.  This  is  because  that  plastic  straining  is  mainly 
concentrated  in  the  ITZs  and  all  aggregates  are  surrounded  by  an  ITZ,  so  that  the  negative  effect  of  the 
volume  fraction  of  ITZ  on  the  average  reaction  stress  level  is  more  significant  than  the  positive  effect 
of  the  volume  fraction  of  aggregates.  However,  it  should  be  noted  that  an  impact  tensile  velocity  V  = 
600  mm/s  and  a  viscosity  parameter  rj  =  400  s~^  are  used  for  the  three  computations.  From  the  former 
results  in  this  section  for  the  influence  of  impact  velocity  and  viscosity  parameter,  it  can  be  expected 
that  when  the  impact  velocity  and  viscosity  parameter  are  decreased,  because  of  a  less  intense  plastic 
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straining  in  the  ITZs,  the  peak  value  of  the  average  reaction  stress  along  the  loading  edge  may  increase 
with  the  increase  of  the  volume  fraction  of  aggregates. 

In  an  attempt  to  simulate  the  fracture  of  the  aggregates  in  the  microstructure  under  impact  tensile 
loading,  another  two  computations  are  conducted.  Microstructure  one  with  mesh  lA  is  used  and  the 
strength  of  the  aggregates  is  reduced  to  the  same  as  the  matrix  material,  i.e.  5.25  MPa.  The  development 
of  plastic  strain  meature  k  is  shown  in  Fig.  5.15  for  the  impact  velocity  V  =  400  mm/s,  600  mm/s  and 
the  viscosity  parameter  rj  =  200  s“F  The  stroboscopic  plots  of  k  along  the  vertical  line  across  the  center 
of  the  aggregate  at  the  lower  left  corner  (see  Fig.  5.1a)  are  given  in  Fig.  5.16.  These  results  show  that 
although  the  plastic  straining  is  still  mainly  concentrated  in  the  ITZs,  a  shear  band  initiated  from  one  of 
the  fixed  corner  does  propagate  through  the  aggregates  and  across  the  whole  microstructure  at  an  impact 
velocity  V  =  600  mm/s.  When  the  impact  velocity  is  reduced  to  V  =  400  mm/s,  this  phenomenon  does 
not  occur.  The  reaction  stress  history  for  the  two  computations  is  shown  in  Fig.  5.16c.  Comparing  with 
the  results  in  Fig.  5.9a,  it  can  be  seen  that  the  peak  value  of  the  average  reaction  stress  does  not  change 
significantly  with  the  decrease  of  the  matrix  strength. 


5.3  Duvaut-Lions  viscoplastic  model 

The  Duvaut-Lions  viscoplastic  model  combined  with  a  Rankine  yield  criterion,  an  associative  flow  rule, 
an  exponential  softening  law  and  a  strain  hardening/softening  hypothesis  is  used  in  this  section.  The 
material  parameters  are  shown  in  Table  5.3  and  different  impact  velocities  V  and  relaxation  times  r 


Table  5.3:  Material  parameters  for  Duvaut-Lions  viscoplastic  model 


Matrix 

Aggregate 

ITZ 

E  3.5e-l-4  MPa 

u  0.3 

CTo  5.25  MPa 

CToo  0.0 

P  0.67e-f3 

p  2.0e-9  Ns^ /mrrP 

E  5.5e-|-4  MPa 

u  0.3 

CTo  41.25  MPa 

CToo  0.0 

P  l.Oe-1-3 

p  2.0e-9  Ns^/mnP 

E  3.5e-f-4  MPa 

V  0.3 

do  3.5  MPa 

dco  0.0 

P  l.Oe+3 

p  2.0e-9  Ns^/mrrP 

Figure  5.17;  Development  of  plastic  strain  measure  k  for  Duvaut-Lions  viscoplastic  model. 
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(a)  T  =  0.5  fj-s 


(b)  T  =  1.0  IJ.S 


(c)  Different  relaxation  times 


Figure  5.18:  Influence  of  the  impact  velocity  (a)  and  (b),  and  the  relaxation  time  (c)  on  the  reaction 
stress  history  along  the  loading  edge  for  Duvaut-Lions  viscoplastic  model. 


(a)  V  =  500  mm/s,  t  =  0.5  fis  (b)  V  =  600  mm/s,  t  =  0.2  /is  (c)  V  =  600  mm/s,  t  =  1.0  /is 


Figure  5.19;  Stroboscopic  development  of  plastic  strain  measure  k  along  the  central  vertical  cross  section 
with  different  impact  velocities  and  relaxation  times. 


(a)  Profiles  of  k  along  the  central  vertical 
cross  section  at  time  t  =  8.0  /is 


(b)  History  of  average  reaction  stress  com¬ 
ponent  CTyy  along  the  loading  edge 


Figure  5.20;  Mesh  objectivity  for  Duvaut-Lions  viscoplastic  model  with  V  =  600  mm/s  and  r  =  1.0  /is. 


are  used  in  the  simulations.  The  microstructure  one  is  analyzed  first.  The  development  of  viscoplastic 
strain  measure  k  with  V  =  600  mm/s  and  t  =  0.5  /is  is  shown  in  Fig.  5.17a  and  b.  Comparing  with  the 
localization  patterns  for  the  Perzyna  model  (cf.  Fig.  5.7)  it  can  be  seen  that  the  highest  plastic  strain  level 
still  occurs  in  the  ITZs  but  a  really  diffuse  mode  of  failure  occurs  in  the  matrix  and  there  are  no  shear 
bands  formed  in  the  matrix.  The  plastic  strain  level  in  the  matrix  is  comparable  to  that  in  the  ITZs. 
The  matrix  material  failed  after  the  reflection  of  the  loading  waves  on  the  fixed  edge,  the  plastic  strain 
increases  consistently  in  the  course  of  time.  These  characteristics  are  due  to  the  choice  of  a  Rankine  yield 
criterion  where  the  maximum  principal  tensile  stress  dominates  the  solution  after  yielding. 

When  the  impact  velocity  is  increased  or  the  relaxation  time  is  decreased  the  plastic  strain  level  is 
higher  (see  Fig.  5.17),  whereas  the  plastic  strain  along  the  fixed  edge  increases  most  significantly.  When 
relaxation  time  r  is  reduced  to  0.2  /is,  the  localization  zone  along  the  fixed  edge  becomes  the  dominant  one 
in  the  matrix  with  the  amplitude  the  same  as  that  in  the  ITZs.  This  can  be  observed  in  the  stroboscopic 
plot  for  plastic  strain  measure  k  in  Fig.  5.19b.  This  mode-I  failure  of  the  microstructure  is  due  to  the 
use  of  the  Rankine  yield  criterion  and  the  absence  of  the  aggregates  along  the  fixed  edge. 


(a)  Structure  one 


(b)  Structure  two 


(c)  Structure  three 


Figure  5.21:  The  deformation  pattern  at  time  t  =  8.0  /ts  for  mesh  lA  (a),  structure  two  (b)  and  structure 
three  (c)  with  V  =  600  mm/s  and  t  =  0.2  /ts. 


Figure  5.22:  Development  of  plastic  strain  measure  k  for  structure  two  (a)-(c)  and  structure  three  (d)-(f) 
with  V  =  600  mm/s  and  t  =  0.2  /xs. 


The  influence  of  impact  velocity  and  relaxation  time  on  the  reaction  force  along  the  loading  edge  is 
presented  in  Fig.  5.18.  It  is  observed  that  the  peak  value  of  average  reaction  stress  component  ayy  is 
not  affected  by  the  impact  velocity  when  a  relaxation  time  t  =  0.5  /xs  is  used,  the  reason  is  that  the 
decrease  of  the  backbone  stress  level  due  to  softening  material  law  neutralizes  the  increase  of  stress  due 
to  the  increase  of  loading  rate.  When  the  relaxation  time  is  increased  to  r  =  1.0  /xs,  the  peak  value  of  the 
average  reaction  stress  increases  with  the  increase  of  impact  velocity.  Combining  the  results  in  sections 
5.1  and  5.2,  it  can  be  concluded  that  the  viscous  effect  plays  a  dominant  role  in  simulation  of  the  positive 
loading  rate  effect  of  the  concrete  microstructure.  The  viscous  effect  can  be  seen  for  different  values  of 
the  relaxation  time  in  Fig.  5.18c,  where  the  peak  stress  level  increases  significantly  when  the  relaxation 
time  is  doubled. 

The  stroboscopic  plots  of  the  development  of  plastic  strain  measure  k  along  the  central  vertical  cross 
section  for  three  numerical  examples  are  given  in  Fig.  5.19.  The  characteristics  of  the  development  of 
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(a)  Structure  two  (b)  Structure  three  (c)  Reaction  stress 

Figure  5.23:  Stroboscopic  development  of  plastic  strain  measure  k  along  the  central  vertical  cross  section 
for  structure  two  (a)  and  structure  three  (b);  The  average  reaction  stress  component  ayy  along  the  loading 
edge  for  the  three  structures. 


(a)  Plastic  strain  measure  k  (b)  Stress  component  ayy 


Figure  5.24:  Spurious  hardening  effect  with  large  plastic  strain  for  V  =  600  mm/s  and  t  =  0.2  fis  at  time 
t  =  5.5  fis. 

plastic  strain  are  exhibited. 

The  mesh  objectivity  is  shown  in  Fig.  5.20  with  the  plots  for  the  profiles  of  plastic  strain  measure  k 
along  the  central  vertical  cross  section  and  the  average  reaction  stress  component  ayy  along  the  loading 
edge  versus  time.  A  finer  mesh  is  needed  in  the  interface  transition  zones. 

The  deformation  of  the  microstructure  amplified  with  a  factor  of  1200  for  V  =  600  mm/s,  t  =  0.2  fis 
and  at  time  8.0  /xs  is  shown  in  Fig.  5.21a.  The  deformation  is  localized  in  the  ITZs  and  the  matrix  near 
the  fixed  edge. 

For  the  microstructures  two  and  three  the  impact  tensile  velocity  V  =  600  mm/s  and  the  relaxation 
time  r  =  0.2  ^s  are  utilized.  The  development  of  plastic  strain  measure  k  is  shown  in  Fig.  5.22.  Because 
of  the  mode-I  localization  pattern  along  the  fixed  edge,  the  additional  aggregates  in  microstructure  two 
have  no  significant  influence  on  this  pattern  and  the  plastic  strain  level.  There  are  more  ITZs  along  the 
fixed  and  loading  edges  for  structure  three,  the  plastic  strain  in  these  ITZs  have  high  values.  However, 
the  plastic  strain  level  in  the  matrix  near  the  fixed  edge  decreases. 

The  stroboscopic  development  of  k  along  the  central  vertical  cross  section  for  these  examples  is  pre¬ 
sented  in  Figs.  5.23a  and  b.  The  plot  in  Fig.  5.23c  compares  the  average  reaction  stress  component  ayy 
along  the  loading  edge  for  the  three  structures  with  mesh  lA  for  structure  one.  We  can  observe  that  the 
average  stress  level  is  nearly  the  same  in  the  beginning.  Then  it  decreases  gradually.  This  is  due  to  the 
reflection  of  loading  waves  on  the  ITZs  with  an  intense  plastic  straining  in  these  three  computations.  As 
the  fraction  of  the  ITZ  increases  (i.e.  from  structure  one  to  structure  three),  the  average  reaction  stress 
level  decreases  during  this  process.  However,  when  the  reflected  waves  on  the  fixed  edge  reach  the  loading 
edge,  the  average  reaction  stresses  for  the  three  structures  increase  to  nearly  the  same  level.  It  indicates 
that  the  negative  effect  of  the  volume  fraction  of  ITZ  on  the  peak  value  of  reaction  stress  neutralizes  the 
positive  effect  of  the  volume  fraction  of  aggregates.  As  for  the  Perzyna  model  in  section  5.2,  it  is  expected 
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(d)  t  =  2.0  fis  (e)  t  =  5.5  fis  (f)  t  =  8.0  /is 

Figure  5.25:  Development  of  plastic  strain  measure  k  in  microstructure  one  with  mesh  lA  for  V  =  400 
mm/s  (a)-(c),  V  =  600  mm/s  (d)-(f)  and  t  =  0.5  fis. 


Figure  5.26:  Stroboscopic  development  of  plastic  strain  measure  k  along  the  vertical  line  across  the  center 
of  the  aggregate  at  the  lower  left  corner  of  microstructure  one  with  r  =  0.5  fis]  The  average  reaction 
stress  component  Cyy  along  the  loading  edge  for  the  two  computations. 


when  the  impact  velocity  is  decreased  or  the  relaxation  time  is  increased  the  peak  value  of  the  average 
reaction  stress  along  the  loading  edge  will  increase  with  the  increase  of  the  volume  fraction  of  aggregates. 
The  deformation  patterns  for  these  two  structures  at  time  t  =  8.0  fis  can  be  seen  in  Figs.  5.21b  and  c. 

It  should  be  mentioned  that  for  simulations  using  the  viscoplastic  models,  we  can  observe  a  physically 
unrealistic  strain  rate  hardening  effect  when  the  plastic  strain  in  some  ITZs  is  very  large.  This  phe¬ 
nomenon  is  especially  significant  for  the  ITZs  located  on  the  fixed  or  loading  edge.  The  plots  in  Fig.  5.24 
provide  an  example,  where  the  arrow  A  points  to  the  ITZ  which  achieves  a  high  plastic  strain  and  stress 
level  at  the  same  time. 

Two  computations  are  carried  out  to  simulate  the  fracture  behavior  of  aggregates.  The  microstructure 
one  with  mesh  lA  is  used  and  the  strength  of  the  aggregates  is  reduced  to  the  same  as  the  matrix  material, 
i.e.  5.25  MPa.  The  impact  velocities  V  =  400  mm/s,  600  mm/s  and  relaxation  time  r  =  0.5  ^s  are  adopted. 
A  varied  smaller  time  step  is  used  for  the  computation  with  V  =  600  mm/s.  The  development  of  plastic 
strain  measure  k  is  shown  in  Fig.  5.25.  The  stroboscopic  plots  of  k  along  the  vertical  line  across  the  center 
of  the  aggregate  at  the  lower  left  corner  (see  Fig.  5.1a)  are  given  in  Fig.  5.26.  These  results  show  that 
aggregate  fracture  occurs  after  the  loading  waves  have  reflected  on  the  fixed  edge.  The  plastic  straining 
in  the  ITZs  is  more  intense  than  that  in  the  matrix  and  the  aggregates.  But  in  the  latter  two  phases, 
the  trends  of  plastic  straining  are  very  consistent.  Compared  with  the  computations  using  the  Perzyna 
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model,  the  aggregate  fracture  still  occurs  when  the  impact  velocity  is  reduced  to  400  mm/s.  Due  to  the 
aggregate  fracture,  the  peak  value  of  the  average  reaction  stress  decreases  when  the  impact  velocity  is 
increased  from  400  mm/s  to  600  mm/s  (cf.  Fig.  5.18a). 


5.4  Constant  gradient  damage  model 

The  second-order  implicit  constant  gradient  damage  model  described  in  section  3.2.1  is  used  in  this 
section.  A  linear  interpolation  for  the  nonlocal  equivalent  strain  is  adopted  in  addition  to  the  linear 
interpolation  for  the  displacement  field.  The  natural  boundary  condition  (Eq.  3.48)  is  enforced  on  four 
edges  of  the  microstructure  by  using  a  penalty  factor  technique.  An  additional  term 

a  f  (5(n^V£<.g)(n^Vee,)dr  (5.1) 

Jan 

where  a  is  the  penalty  factor,  is  augmented  to  the  weak  form  of  the  diffusion  equation  (Eq.  3.59).  The 
exponential  softening  damage  evolution  law  and  the  modified  Von  Mises  definition  of  equivalent  strain 
are  utilized.  The  material  parameters  are  shown  in  Table  5.4. 


Table  5.4:  Material  parameters  for  implicit  constant  gradient  damage  model 


Matrix 

Aggregate 

ITZ 

E 

3.5e+4  MPa 

E 

5.5e-l-4  MPa 

E 

3.5e-l-4  MPa 

V 

0.3 

1/ 

0.3 

V 

0.3 

Ki 

1.5e-4 

7.5e-4 

l.Oe-4 

Q 

0.96 

a 

0.96 

a 

0.96 

fi 

500.0 

P 

500.0 

p 

500.0 

p 

2.0e-9  Ns^/mm^ 

P 

2.0e-9  Ns^/mm^ 

p 

2.0e-9  Ns^/mrrP 

k 

1.0 

k 

1.0 

k 

1.0 

(f)  u),  t  =  8.0  fis 


Figure  5.27:  Development  of  equivalent  strain  (a)-(c)  and  damage  (d)-(f)  with  V  =  400  mm/s  and  c  = 
0.05  mm^  for  constant  gradient  damage  model  and  mesh  lA. 


The  development  of  equivalent  strain  e^g  and  damage  uj  with  impact  velocity  V  =  400  mm/s  and 
gradient  parameter  c  =  0.05  mm^  for  structure  one  mesh  lA  is  presented  in  Fig.  5.27.  We  observe  that 
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(a)  Different  impact  velocities  (b)  Different  impact  velocities  (c)  Different  gradient  parameters 

Figure  5.30:  Influence  of  impact  velocity  and  gradient  paramter  on  the  reaction  stress  history  along  the 
loading  edge  for  constant  gradient  damage  model  with  mesh  lA. 


y  (mm)  y  (mm)  y  (mm) 

(a)  V  =  400  mm/s,  c  =  0.05  mm^  (b)  V  =  500  mra/s,  c  =  0.05  mm^  (c)  V  =  400  mm/s,  c  =  0.02  mm^ 

Figure  5.31:  Stroboscopic  development  of  equivalent  strain  Ceq  along  the  central  vertical  cross  section  of 
the  structure  one  mesh  lA. 


y  ("""j  t  (•) 

(a)  Profile  of  te,  along  the  central  vertical  (b)  Average  reaction  stress  component  cryy 

cross  section  at  time  t  =  8.0  /xs  history  along  the  loading  edge 

Figure  5.32:  Mesh  objectivity  of  constant  gradient  damage  model  for  structure  one  with  V  =  400  mm/s 
and  c  =  0.05  mm^. 

enhanced  models  is  irrespective  of  the  size  of  imperfection.  Note  that  the  mode-II  failure  pattern  is  due 
to  the  fact  that  when  k  =  1,  the  definition  of  equivalent  strain  is  a  standard  Von  Mises  model. 

The  impact  velocity  increases  to  V  =  500  mm/s  in  a  simulation  presented  in  Fig.  5.28  with  the  same 
material  parameters  as  that  in  Fig.  5.27.  Before  the  refiection  on  the  fixed  edge,  the  loading  waves  cause 
much  more  damage  in  the  ITZs.  After  reflection  a  similar  failure  process  as  that  in  Fig.  5.27  is  observed, 
while  the  strain  level  and  damage  level  are  much  higher. 

When  the  impact  velocity  is  still  taken  as  V  =  400  mm/s  but  the  gradient  parameter  is  reduced  to 
0.02mm^,  there  are  more  damage  and  deformation  occurring  in  the  microstructure  as  that  for  an  increased 
impact  velocity.  The  evolution  of  equivalent  strain  and  damage  is  given  in  Fig.  5.29.  The  failure  pattern 
remains  unchanged  but  the  width  of  localization  zones  becomes  narrower. 


(a)  Development  of  €cq  for  structure  three 


(b)  Reaction  stress  history  for  the  three 
structures 

Figure  5.35:  The  stroboscopic  development  of  equivalent  strain  Cgg  along  the  central  vertical  cross  section 
of  the  structure  three  (a)  and  the  average  reaction  stress  component  cTyy  along  the  loading  edge  versus 
time  for  the  three  structures  (b). 


(a)  Structure  one,  mesh  lA  (b)  Structure  two  (c)  Structure  three 

Figure  5.36:  Deformation  patterns  for  the  three  structures  with  V  =  400  mm/s,  c  =  0.05  mm^  and  at 
time  t  =  8.0  /iS. 


viscoplastic  strain  rate.  It  can  be  concluded  that  at  relatively  high  impact  tensile  velocity,  a  positive 
loading  rate  dependence  can  hardly  be  obtained  with  this  gradient  dependent  model. 

The  development  of  equivalent  strain  along  the  central  vertical  cross  section  of  microstructure  one  with 
mesh  lA  is  presented  by  using  stroboscopic  plots  in  Fig.  5.31.  The  above  analyzed  characteristics  of  the 
numerical  results  can  also  be  observed  in  this  figure. 

The  mesh  objectivity  of  the  numerical  solution  with  the  use  of  a  constant  gradient  damage  model 
is  checked  by  comparing  the  results  obtained  for  the  two  meshes  of  structure  one.  The  distribution  of 
equivalent  strain  along  the  central  vertical  cross  section  of  the  microstructure  at  time  t  =  8.0  /iS  and 
the  average  reaction  stress  component  ayy  along  the  loading  edge  versus  time  are  ploted  in  Fig.  5.32.  It 
shows  that  the  solution  is  independent  of  the  finite  element  discretization. 

A  similar  numerical  simulation  is  carried  out  on  microstructure  two  and  microstructure  three  with 
impact  velocity  V  =  400  mm/s  and  gradient  parameter  c  =  0.05  mm^.  The  evolution  of  equivalent 
strain  e^g  and  damage  a*  for  structure  two  is  shown  in  Fig.  5.33.  It  shows  that  the  additional  aggregates 
in  the  microstructure  deviate  the  propagation  of  the  dominant  shear  band.  The  strain  level  does  not 
change  significantly  compared  with  Fig.  5.27,  but  there  are  more  steep  strain  gradient  variations  in  the 
microstructure  due  to  the  larger  volume  of  ITZ.  The  aggregate  volume  is  further  increased  in  structure 
three,  so  is  the  volume  of  ITZ.  The  numerical  results  are  shown  in  Fig.  5.34.  Because  of  a  very  high 
strength,  the  aggregates  do  not  fail  during  loading.  They  obstruct  and  deviate  the  propagation  of  local¬ 
ization  zones.  Hence,  a  diffuse  mode  of  failure  is  formed  for  this  structure.  The  stroboscopic  development 
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of  equivalent  strain  is  given  in  Fig.  5.35a  for  structure  three. 

The  comparison  of  the  average  reaction  stress  history  along  the  loading  edge  for  the  three  structures 
with  the  same  material  parameters  and  impact  velocity  is  shown  in  Fig.  5.35b.  The  limit  stress  level 
increases  with  the  volume  of  aggregates.  Comparing  with  the  results  in  sections  5.2  and  5.3,  it  can  be 
concluded  that  the  influence  of  a  higher  volume  of  aggregates  is  more  significant  when  the  damage  pattern 
is  more  diffuse  and  this  strengthening  effect  is  greatly  reduced  when  the  damage  is  mainly  concentrated 
in  the  ITZ,  unless  a  viscous  effect  is  incorporated  in  the  model.  The  deformation  patterns  for  the  three 
structures  at  time  t  =  8.0  fis  are  shown  in  Fig.  5.36. 
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Chapter  6 
Discussion 


In  this  report  different  continuum  constitutive  models  are  used  to  investigate  the  dynamic  response  of 
concrete  at  meso-level  under  impact  tensile  loading. 


For  the  one-dimensional  bar  problem,  mesh  dependence  of  the  numerical  results  as  one  of  the  con¬ 
sequences  of  the  ill-posedness  of  the  initial  value  problem,  is  illustrated  with  a  standard  plastic  model 
combined  with  a  linear  softening  law  and  a  local  damage  model.  Incorporation  of  higher-order  temporal 
derivatives  or  higher-order  spatial  derivatives  into  the  continuum  model  can  regularize  the  governing  field 
equations  and  preserve  the  well-posedness  of  the  initial  value  problem  upon  the  onset  of  strain  localiza¬ 
tion.  The  Perzyna  viscoplastic  model,  the  Duvaut-Lions  viscoplastic  model  and  the  gradient  enhanced 
damage  model  are  used  to  exhibit  the  regularization  effect.  In  addition  to  a  constant  gradient  damage 
model,  a  strain-based  transient-gradient  damage  model  is  also  used  to  illustrate  that  making  the  gradient 
parameter  a  variable  can  avoid  suprious  spreading  of  damage  at  the  final  stage  in  a  failure  analysis. 


Three  block  square  microstructures  are  used  to  conduct  two-dimensional  analyses.  Three  phases  are 
considered  in  the  microstructure,  namely  matrix,  aggregate  and  interfacial  transition  zone.  Physically 
meaningful  result  can  not  be  obtained  by  using  a  standard  plastic  model  because  the  numerical  results 
show  an  obvious  mesh  dependence.  When  a  Perzyna  model  combined  with  a  Von  Mises  yield  criterion 
is  used,  two  types  of  localization  zones  can  be  distinguished.  One  type  includes  those  initiated  from  the 
ITZ  and  which  are  controlled  by  the  width  of  the  ITZ,  while  the  other  type  consists  of  those  formed  in 
the  matrix  and  are  dominated  by  the  internal  length  scale  of  the  material.  The  plastic  strain  level  in 
the  ITZs  is  much  higher  than  in  the  matrix.  The  localization  zones  formed  in  the  matrix  can  propagate 
but  can  not  develop  into  dominant  shear  bands  due  to  the  constraints  imposed  by  the  aggregates.  The 
plastic  deformation  increases  with  the  increase  of  impact  velocity  and  decrease  of  a  viscous  effect.  A  finer 
discretization  is  needed  in  the  ITZs  in  order  to  capture  large  strain  gradients.  The  confining  effect  of 
aggregates  on  the  localization  zones  in  the  matrix  is  further  illustrated  by  the  simulations  on  structure 
two  and  three.  In  the  computations,  a  positive  loading  rate  dependence  is  obtained  which  is  due  to  the 
viscous  effect  incorporated  in  the  constitutive  relation.  The  peak  value  of  average  reaction  stress  along 
the  loading  edge  decreases  with  the  increase  of  volume  fraction  of  the  aggregates.  However,  it  is  expected 
that  the  enforcing  effect  of  aggregates  can  be  simulated  when  the  viscous  effect  is  increased  in  the  material 
model.  Additional  simulations  with  a  reduced  strength  for  the  aggregates  show  that  the  mode-II  fracture 
of  aggregates  can  be  simulated  with  the  Perzyna  model. 

When  a  Duvaut-Lions  viscoplastic  model  combined  with  a  Rankine  yield  criterion  is  used,  a  localization 
pattern  which  is  different  from  that  for  the  Perzyna  model  is  obtained.  The  plastic  strain  level  in  the 
matrix  is  comparable  with  that  in  the  ITZs,  and  it  increases  consistently  with  loading.  Hence,  a  more 
diffuse  failure  mode  can  be  formed.  However  when  the  impact  velocity  is  increased  or  the  relaxation  time 
is  reduced,  a  dominant  localization  zone  of  mode-I  type  emerges  along  the  fixed  edge.  Wave  reflections 
on  the  aggregates  can  cause  earlier  and  more  intense  plastic  straining  in  the  matrix  and  the  ITZs.  At  the 
same  time,  the  increase  of  the  aggregate  volume  can  increase  the  average  modulus  and  the  load  carrying 
capacity  of  the  structure.  These  two  effects  can  be  observed  in  the  numerical  results.  The  peak  value  of 
average  reaction  stress  exhibits  no  loading  rate  dependence  with  a  specific  choice  of  the  relaxation  time. 
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However,  when  the  relaxation  time  is  doubled  a  positive  loading  rate  dependence  is  obtained.  The  same 
trend  can  be  observed  for  the  simulation  of  the  enforcing  effect  of  aggregates.  When  the  strength  of  the 
aggregates  is  reduced,  a  mode-I  fracture  of  aggregates  can  be  observed  in  the  numerical  results. 

The  third  model  adopted  is  the  second-order  implicit  constant  gradient  damage  model.  Utilizing  the 
modified  Von  Mises  definition  of  equivalent  strain,  a  dominant  shear  band  can  be  formed  in  microstructure 
one.  The  intensity  of  strain  localization  increases  with  the  increase  of  impact  velocity  and  the  decrease  of 
gradient  parameter.  The  band  width  decreases  when  the  gradient  parameter  is  reduced.  The  additional 
aggregates  in  microstructure  two  change  the  propagation  of  the  dominant  shear  band,  whereas  due  to 
the  high  volume  of  aggregates  in  microstructure  three  the  localization  pattern  is  more  diffuse.  In  the 
computations,  the  peak  value  of  the  average  reaction  stress  decreases  with  the  increase  of  the  impact 
velocity.  A  positive  loading  rate  dependence  can  hardly  be  obtained  at  relatively  high  impact  velocity 
due  to  the  absence  of  a  viscous  effect.  When  the  volume  fraction  of  aggregates  is  increased  the  maximum 
reaction  stress  level  along  the  loading  edge  increases. 
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